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Scope 


This is the Annual Report prepared by the Smithsonian Astrophysical Observatory 
for the second year research activity on NASA Marshall Space Flight Center Grant 
NAG8-1046 with Charles Rupp as Technical Monitor. This report covers the period 
from 22 February 1995 through 21 February 1996. 
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Summary 


This report deals with the following topics which are all related to the orbital 
injection of the SEDSAT satellite: 

1 ) Dynamics and Stability of Tether Oscillations after the First Cut 

The dynamics of the tether after the first cut (i.e., without the Shuttle attached to it) is 
investigated. The tether oscillations with the free end are analyzed in order to assess the 
stability of the rectilinear configuration in between the two tether cuts. 

2) Analysis of Unstable Modes 

The unstable modes that appear for high libration angles are further investigated in 
order to determine their occurrences and the possible transition from bound librations to 
rotations. 

3) Orbital Release Strategies for SEDSAT 

A parametric analysis of the orbital decay rate of the SEDSAT satellite after the two 
tether cuts has been carried out as a function of the following free parameters: libration 
amplitude at the end of deployment, deviation angle from LV at the first cut, and orbital 
anomaly at the second cut. The values of these parameters that provide a minimum 
orbital decay rate of the satellite (after the two cuts) have been computed. 

4) Dynamics and Control of SEDSAT 

The deployment control law has been modified to cope with the new ejection 
velocity of the satellite from the Shuttle cargo bay. New reference profiles have been 
derived as well as new control parameters. Timing errors at the satellite release as a 
function of the variations of the initial conditions and the tension model parameters have 
been estimated for the modified control law. 
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INTRODUCTION 


Some aspects of the dynamics of SEDSAT orbit injection are studied. After being deployed from the space 
Shuttle flying in a circular orbit, the SEDSAT satellite and tether will librate zenithwards as a gravity gradient 
pendulum [1], At a suitable moment, the tether will be cut at the orbiter end, leaving the tethered system (TS) 
to enter a higher energy elliptical orbit, safely away from the shuttle. If the TS is straight at cut time, it will 
continue librating about its center of mass (COM) keeping its rectilinear shape, just as a rigid body. A second 
cut is done, leaving the satellite and tether to continue separate orbits. The libration can be used to transfer 
angular momentum from the tether to the satellite, which will reach a higher orbit at the expense of the tether. 
The mass of the tether is not negligible. The four stages — deployment, libration, TS orbit and satellite orbit — 
are to be optimized to obtain the longest possible satellite life. 

The COM orbit sensitivity to deployment parameters and cut time, studied by Jesus Pelaez in Ref. [2], 
will be used along this research. 

The main topics to be studied are the stability of the rectilinear configuration and the optimization of 
operational life. 

It is known that a straight shape is a solution of the equations of motion of the TS [3, p. 39]. This is 
precisely the movement desired in order to obtain the maximum transfer of angular momentum from the tether 
to the satellite. However, we must study whether this movement is stable, since it will not be possible to get 
the perfect initial conditions at cut time. On the contrary, the tension released at cutting will result in tether 
shrinking, producing speed and position deviations from the rectilinear shape. We need to know how the 
perturbations of the rectilinear shape will propagate during the few orbits the injection will take. The main 
danger comes from the free end. The stability of two-mass tethered satellites has been thoroughly studied: the 
tension produced by rotation and the two end masses has a stabilizing effect. But in this case one tether end is 
free [4, 5]. The instability is more likely to appear when the TS swings backwards in the libration, at which 
time the tension falls to its minimum. Besides, most stability studies have been done for circular orbits, while 
here the elliptic orbit will induce a forcing term from the inertia forces. The oscillation modes and frequencies 
will be sought by different methods and the evolution of the libration itself will be investigated. 

Another point of concern is selecting the right set of parameters to obtain the maximum satellite useful 
life. Since the shuttle orbit is fixed, we have only three free parameters: libration amplitude, libration angle at 
first cut time, and second cut time. The maximum operational life corresponds to maximum orbit energy and 
minimum eccentricity. This corresponds to minimum orbital decay rate. We will compute the orbital decay rate 
as a function of the three free parameters. 

These two aspects are closely related. It is obvious that the wider the libration amplitude, the greater the 
speed increase obtained at cut time, and the higher the orbit attained. However, these large librations bring us 
closer to the instability region. We need to know how far we can go without a catastrophic propagation of 
perturbations. 
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Stability of Pendular Motion 


Elliptic Orbit of the Center of Mass 

Within the degree of approximation considered in this study, the center of mass of the TS will follow an 
unperturbed elliptic orbit. The Deployment, first cut, and second cut will take just a few orbits, so that no 
perturbations need to be considered: only spherical earth gravitation and inertia forces will be considered. The 
parameters of the COM orbit are known as a function of the deployment and cut data [2]: 
a=R( 1 +2( 1 +u)e+. . . ) ; 
e=2ue+...; 

<#=oob(l+3(l+u)e+...); 

Where u=3+2V3sin0 M and e=^/2Ro. 8 M is the TS libration amplitude; Ro is the shuttle’s circular orbit 
radius; ^ is the distance from the free end to the center of mass of the TS with the tether straight, while Lf is 
the final deployed tether length. We will use here the same notation as [2] for u and e, but a different one will 
be used later on for the study of decay optimization in section 3. 

If the TS is straight at cut time, it will keep its straight shape and orbit as a rigid body [3], The terms 
coupling COM motion with rotation about the COM are smaller than the approximation considered here, so 
that the gravidyne effect studied by Beletski is completely negligible for the size of the TS considered [3, pp. 
67-68], We can thus study separately the motion of the COM and the TS motion about its COM. 


Pendular Motion Equations 



Figure 1 Reference frame for 
theTS. 


The tether is considered as perfectly flexible and inextensible. A 
reference frame Gxyz is chosen as follows: origin in the tethered system 
center of mass, Gx aligned with the local vertical towards zenith, Gz 
aligned with the orbital speed of the center of mass, and Gy completing a 
right-handed frame. The equations of the TS motion are: 



J 0 'pRds+m A R A =0 


Where t is time and s the curvilinear abscissa, with origin at the end 
mass; p is the linear density of the tether, U its total length, T(s,t) its 
tension, f the force per unit mass, R(s,t) the position vector, and A the 
end mass. The second equation corresponds to a perfectly inextensible 
tether, and the third to taking the COM as origin of the reference frame. 
We shift to non-dimensional variables using the average orbital rate and 
tether length: 


x = n(t-t p ), s = L f ct, R = L f r, T = pn 2 L J f F; 


where n = 



Obtaining: 


r = (Fr')' + 


L f n 2 


M=i 

£r(a, i)da 


+- — -r(0, t) = 0 
m. 


( 2 ) 


With boundary conditions: 
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a = 0 


dF pL f 

|L(0 > x) = i^-!-F(0,T) 
da m A 

F-r-r-(0,x) = 0 


da 


a = l F = 0 

Where 0=0 is the end mass and 0=1 the free end. The perturbing force is: 

. r — 3xi 


f =-— 


+- 


ear 


xcosv- 2zsin v 


z 

0 

• + 2©< 

0 

zcosv + 2xsinv 


-X 


Where v is the true anomaly of the COM orbit; r=(x,y,z); K=l+ecosv; co = v. The 1/Lf factor in f can be 
dropped if x,y,z are non-dimensional lengths. Only terms of order e and ee are kept; e 2 , ee 2 , and e 2 e have been 
dropped. 

If, hypothetically, the tether is assumed to keep a rectilinear shape, it will be sufficient to give the position of 
one of the masses in spherical coordinates [3, p. 378ff]. The equations of motion of a fully deployed massive 
tether, considered as a rigid body, and without any perturbations nor thrust (dumbbell model), are: 

0 + (b y - 2<p tan cp^0 + co y ) +■ -^sin0 cos 0 = 0 

k rz 


<P + 
T = 


k g 

(^“yf-irl 1 - 300520 ) 


sin (p cos (p =0 


(3) 


P(s) 


P(s) = m A (s G -s A )+£p(s)(s G -s)ds 


Where 0 and cp are the spherical coordinates of A. Note that the reference frame used here is different from 
the one in [3]. These expressions are valid only while T>0. The tension is a function of t and s, which, for 
constant p and only one mass, is: 

s 2 0 2 L 2 

P(s) = -p y + ps c S + m A s c = -p + pbL 2 f a + m A b L, 

s 1 

Where b is the non-dimensional curvilinear abscissa of the COM: b = — = — 


1 + 


pL f 


We shift to the non-dimensional arc length v, which is the same as 0 but for taking now G as origin, v is 
positive towards A and is defined as follows: v=b- 0 . Consequently, 




(4) 


If we consider in-plane pendular motion only, (p and its derivatives will be zero, and equations (3) reduce only 
to the 0 equations, as follows: 

Kn 2 0-2co 2 esin v + 3oo 2 sin 0 cos 0 = 0 


2F = 


CO . 

0 + - - 


CD 


Kn 


(l-3cos 2 0) [(b-l) 2 - v 2 j 


(5) 


Where dots are derivatives with regard to non-dimensional time x, which is equivalent to the mean anomaly 
of the COM orbit, co is the instant angular speed of the orbital frame, which corresponds to the orbit of the 
COM. Note that, although we are using non-dimensional time, mis still the derivative of v with regard to real 
time t. We have: 


CD 


_ dv _ nr 

dt v a 3 


(l+ecosv)~ 2 P /, .dm 
3T-; 0)2 = -rV (1 + e cos v); — 

(1 -ef R ° dt 




ecosv) ^ . co^ 

— e sin v = -2e smv — 

72 K 


(l e 2 ) 
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Numerical integration of the equations above requires solving Kepler’s equation at every step. It is therefore 

easier to use the true anomaly as variable: 

• 2esinv/ A \ 3 . A ^ 

0 (0 + 1 ) + — sin 0 cos 0 = 0 

V* ' ' 1C 

( 6 ) 


F = 


00 

2 n 2 


(e + 1 ) (l-3cos 2 0) [(b-l ) 2 -v 2 ] 


Where 9 = 30/ 3v . These equations correspond to the gravity gradient pendulum, plus small damping and 
forcing terms due to inertia forces in an eccentric orbit. The initial conditions for the integration correspond 
to the libration in the circular orbit of the shuttle. An adjustment has to be made to the initial deviation and 
angular speed when the system switches from the circular orbit frame to the new elliptic orbit frame after 
cutting the tether. 


Deviations from the Pendular Motion 

In order to study the deviations from a pendular motion, we will erect a reference frame tied to the assumed 
rigid shape of the TS: origin in the COM, GS, along the tether towards A, Gq parallel to Gy, and G£ normal 
to the tether within the orbital plane. GJ;q£ moves according to the equations above. The coordinates of the 
TS points in the new frame will be: 


(7) 




x = £cos 0 -£sin 0 

ti(ct, x) 

► < 

y = 7i 

C(M. 


z = J;sin 0 + £cos 0 


In the new frame, the derivatives to be used in the TS motion equations will be: 


r = 


"i 


A 




X 

•ii 

■ + 20 < 

0 

► + 0 < 

0 

-e 2 i 

0 

C 
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A. 


X 


(Fr') = 


K' 

Ft|' 

K 


( 8 ) 


ficn 2 = co 2 - 

3^ cos 2 0-^sin0cos0) 

> + eco 2 cos v< 

X 

0 

+ 2eco 2 sin v< 

0 

► + 2mcco< 

'4+0' 

o • 


3(-£cos0sin0 + £sin 2 0) 


X 


A. 


- 4+0 


If we consider the initial conditions of a pendular motion — E,(a,0)=a A -a, q=0, £=0, and first derivatives zero 
except £,’=- 1 — and substitute them into the above equations, we obtain the initial acceleration: 


|(a,0) = -F' + (a A -a)| 
q(a, 0 ) = 0 
C(o,0) = (a A -a)| 


3cd 2 


eco 


2co ‘ 


Kn 


cos 0 + — —cos v + — 0 + 0 


Kn 


(9) 


x 3co 2 . 2eco 2 . ^ 

0 r3cos0sin0 + -smv 


Kn 


Kn 


The initial accelerations are linear in a A -a except for F\ If the tether is initially straight, pendular motion 
equations apply to the TS-body frame and, therefore: 

( « /a 2 \ 

F' = (a A -a) 

Yn _ lm ” n j 

( 10 ) 


Kn 


rar 


3co 2 2eco 2 

0 H -3cos0sin0 —sin v = 0 


Kn 


Kn" 


F’ in (10) is obtained by derivating ( 6 ) with regard to v and then switching back to o. This leads to: 
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^ = ii = C = 0 (11) 

Pendular motion is thus a solution of the equations of movement for the TS when the initial conditions are 
rectilinear. What has to be studied now is the stability of this solution, and the speed at which perturbations 
propagate or die out. 


Small Transverse Oscillations 

Since the equations of the in-plane and out-of-plane deviations £,t|,£ are coupled through the tether continuity 
equation |r'| = 1 and tension only, they can be studied separately. If they are assumed to be small, they can be 

decoupled since the coupling terms are of a higher order. Taking into account the equations for pendular 
motion, we obtain: 


S = (FV) + 

Ti = (Fti')'- 

C = (F?') , + 

%' 2 + T |' 2 


3co' 


' 

^ Kn 
CO 2 


Kn 2 


f 3co 

+ 



v Kn : 

+c 2 = 


2 eco 2 2 co 0 

■cos 0 + —COS V H f-0 

Kir n 


■ 2 n eco 2co0 ■ 2 

-sin 0-1 — cos v H 1-0 


A 


n 


( 12 ) 


Kn" 


-<!+•* 


Plus the condition that G coincides with the COM: 


f rdv + — — r(b) = 0 
L ' PL, ^ ' 


(13) 


Where v=b is the end mass A and v=b-l (negative, since b<l) is the free end. To simplify the notation, the 
coefficients of ^ and £ will be called G(x) and H(t). G is the same function of t appearing in the non- 
dimensional tension F(v,x). If we consider now the deviations from rectilinear shape and pendular motion to 
be small, we have: 

£(v,t) = v+e£,(v,x) 

Ti(v,x) = £n,(v,x) ^e«l (14) 

£(v,t) = ££i(v,t) 

Where we have in fact made a change of variable from a to v=b-a. (’) refers now to derivation with regard to 
v, but signs in (12) remain unchanged due to the double derivation. From the continuity equation we deduce 
that longitudinal perturbations are much smaller than transversal ones: 

l+2^+e^ 2 +e 2 Ti; 2 +£ 2 C; J =l => -2^=e& 2 +r\[ 2 +tf) (15) 

Both £1 and are one order of magnitude smaller than the transverse oscillations. Therefore, we should 
redefine them to show their relative orders of magnitude: 

£(v,x) = v+£ 2 S,(v,x) 

Tl(v,x) = en 1 (v,x) 

C(v,t) = eC(v,t) 

F = F 0 + eF, = G(x)p(v) + eF, (v, x) 

Substituting this into the deviation equations: 


(17) 


(16) 


e 2 ^, =[(F 0 +eF ! )(l + £ 2 ^)] +G(x)(v + £ 2 ^)+2^+0 Je; i 
' (0 2 

Eti, =[(F 0 -t-£Fi)(enO] — r^i 

L Kn 

«£, =[(F„ +£F,Xec;)]' +H(xK, - 2 ^+ej t% 

And keeping only terms of order e: 
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( 18 ) 


o=F;+2^+e)c, 

*i. = ( F JiO 

C.“(f.£) ,+ h(tK, 


While the 0(1) term in the first equation is an already known result: F 0 ' = -vG(x) 
The boundary conditions are: 


|F I '(b) = -^F 1 (b) 


rircb) = o 


m. 


F,(b-1) = 0 


\L"' 


, m A 
dv+- A 


pLf 


rii(b) = 0 


cr(b)=o 

|J>v + ^( b ) = 0 

l J o-i pL, 


(19) 


No condition can be imposed on the free end except that tension be zero. 

As we can see, in the first approximation the transverse in-plane and out-of-plane oscillations are decoupled, 
and longitudinal oscillations are not present. We will continue developing them in parallel, although they 
could be treated separately. 

The first approach to be tried is separation of variables. We assume solutions of the form: 

T 1i( v ,'t) = A(v)Z(T); C,(v.t) = B(v)E(x); F 0 =G(t)P(v) (20) 


A| Z + Z I = GZ(PA') 
Kn 


b(e-h(t)e)=ge(pb') 


z+ 


CO 

K31 2 


( PA 0 


GZ 


=-x. 


(21) 


E-H(t)E (PB') 

GE B " c 


We obtain a Sturm-Liouville problem. We will try to calculate the eigenvalues and eigenfunctions. They will 
be the same for in-plane and out-of-plane oscillations, since the spatial problem is the same for both. 


Analytical Eigenvalues and Eigenfunctions 


The equation to be solved for both transverse oscillations is: 

FA" + F'A' + XA = 0; if we substitute F = j^(b — l) 2 - v 2 J/ 2 , we obtain: 

[(b-1) 2 - v 2 ] A"-2vA' + 2XA = 0 
The boundary conditions are: 

• v=b f£( b)=0 

• v=b-l No condition, since the tip is free. F is already zero there. 

• COM J' A(v)dv + -^-A(b) = 0 

b 1 P^r 

We can substitute = - — — . Besides, the two boundary conditions 

pL r 2b 


( 22 ) 

(23) 

(24) 

turn out to be the same when the 


differential equation is used to calculate the integral; this is expected since the eigenfunctions will have one 
multiplicative constant left free. Besides, since we are studying the movement about the COM, if the initial 
deviations and speed meet (24), it will be fulfilled for all t. We could use either form of the boundary 
condition, or a third one obtained from the differential equation: 

A '( b ) = " A ( b ) (25) 

With the change v=(l-b)x, and taking 2^=n(n+l), we obtain the Legendre Equation: 

[l-x 2 ]A"-2xA' + n(n + l)A = 0 With B.C.: A'|-p-j = 0 (26) 


8 



Solutions will be Legendre functions P n (x) and Q n (x). However, Legendre polynomials of the 2nd kind are not 
bounded in x=-l, and thus cannot be part of the solution. But P n (x) cannot be part of the solution either, since 
they do not meet the boundary conditions: their second derivative becomes zero at the origin, not at x=b/(l-b). 
But b«l, so that the odd-order Legendre polynomials can be a first approximation of the solution. The even- 
order polynomials do not meet the approximate boundary conditions. 

In an unpublished work, Jesus Pelaez and Manuel Ruiz obtained an approximate solution by asymptotic 
expansion of the solution in powers of £ = b / (1 - b) « 0.085 : 


A = A 0 + £ A j H — 
A, = X Q + £A. t 


with: 


A 0 (x) = C 0 P“(x) 
2X 0 — n(n + 1) 


as first approximation, and: 


n(n + l)-2, . 

K =-2 V 7 - , (l + 2n)- 
1 7 tn(n + l) v ' 


1 +- 


1 n 
— + — 

2 2 


(27) 


A, =C,P„°(x)+4^-d 


H 1 + 


n -+- 
2 2 


■Q:W + 2X 1 (Q:(x)j;Q:(y)P:(y)dy + P:(x)j;p;(y) 2 dyj 


as the second. T is the Gamma function. The singularity at the free end due to Q n (x) is controlled because its 
coefficient tends to 0 as x approaches -1. The approximate eigenfunctions (Ao only) are shown in Figure 2. 
Fortunately, there is no need to compute A t for our purposes. With X\ we can compute the time evolution of 
the different eigenmodes. The eigenvalues we will consider are: 



l 

2 

3 

4 

5 

6 

n 

1.0000 

5.3864 

13.3127 

24.7645 

39.7405 

58.2405 


1 

0.8 
0.6 
0.4 
0.2 
0 

- 0.2 
- 0.4 
- 0.6 
- 0.8 
-1 

-1 - 0.8 - 0.6 - 0.4 - 0.2 0 

Non-dimensional Tether Length 


Table 1 


Approximate Eigenfunctions 



Figure 2: Six first approximate Eigenfunctions (A 0 only). 
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Numerical Computation of the Eigenmodes’ Time Equations 


Substituting these values of X in the time equations (21 ), we obtain: 


Z + 


f co 2 ' 

Kn 


Z = 0 


(28) 


E + (W3(x)-H(x))E = 0 

Where H and G are functions of x through v and 6. Dots mean derivation with regard to non-dimensional 
time x. Since 0 is only known as a function of time by integrating (5), there is no way to solve these equations 
analytically. A numerical integration has to be carried out, in parallel with the equations for 0, and this for 
each eigenvalue we want to explore. 

For this, the system being linear, we will take two independent solutions and integrate them for a number of 
orbits to see their evolution. One will have initial deviation and zero initial speed, and the other initial speed 
and zero deviation. We will take the first six eigenvalues, corresponding to n=l,3,5,7,9,l 1 . 

To avoid the difficulty of solving Kepler’s equation at each integration step, we will take v as the independent 
variable. This change implies: 


dx _ dx dv _ . co d 2 x 

dx dv dx n ’ dx 2 
which yields the system of equations: 


_d_ 

dx 


. co 
x— I = x- 


.. co 2esinvco 2 


n 


Kn 


(29) 


0 = 


2esin v 


f© -f- 1) sinGcosO 

\ / K 


g= 2esinvj_ 


i= 2e 2 nv z _ 


{X-l)G + 


M3+- 

K 


3(cos 2 0-sin 2 0) 


(30) 


But now, G and H are functions of v: 

(3cos 2 0 -l) 


G(v ) = (0 + 1) 

H(v) = (0 +l) 2 


(3 sin 2 0-1) 


(31) 


And the initial conditions will be: 

0(0) = 0 E,(0) = 0 E 2 (0) = 1 

0(0) = Q 0 E,(0) = 1 E 2 (0) = 0 


(32) 


Z|(0) = 0 Z 2 (0) = 1 
Z,(0)=1 z 2 (o) = o 

The initial libration angular speed is given by the gravity gradient pendulum speed at the crossing of the local 
vertical, obtained from the energy equation: 0 O = co o V3sin0 M , where Cflb is the circular orbit angular speed 

and 0 M the maximum libration amplitude before cutting. These angle and angular speed are referred to the 
circular orbit frame. After the cut, we use a frame tied to the elliptic orbit of the COM; therefore, speeds have 
to be adjusted. If system 1 is a fixed frame centered in the Earth, 2 the TS-bound frame, 0 the COM orbital 
frame, and 3 the circular orbit frame, the angular speed we seek is CO20, and we can write: 

0(0) = CDj,, = (Ojj + to,, - co 0I ; where 0^3 is 0 O = co 0 a/ 3 sin0 M above, (O31 is the circular angular speed coq, and 
Obi is the angular speed of the elliptical orbit, which we take from [2]: 


co 0 . = v(v = 0) = n 


(l+ecos v) 


(l-e 2 ) 

Adding the terms we obtain: 
d0 




l + 3(l + u)e+.. 


(l +e )' 


(l-e) 3 ' 


= o) 0 [l + (u-3)e+...] (33) 


= — (O) = co o A/3sin0 M (l-2e+...) 


(34) 


with respect to time, and: 
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( 35 ) 


0(0) = — (0) = ^ = V3 sin 6 m fl - 2(V3 sin 0 M + l)e+.. .1 
dv v 0 L v ' J 

with respect to true anomaly, which is the one we will use in the numerical integration. This is done by the 

routine “TRVSA100.C”, which computes the amplitudes of the oscillations for 100 periods of each 
independent solution and eigenvalue (not orbital periods) for values of the initial libration amplitude (0 M ) 
ranging from 0 to 60 degrees. The period of each oscillation is computed as well. 

Remarkable points of the results are: 

• Eigenvalue 0 diverges fast for both in-plane and out-of-plane oscillations. The greater 0 M , the faster it 
grows. This is not as troublesome as it looks, since the corresponding eigenmode, as shown in Figure 2, 
is a straight line. It only reflects the stability of the rigid body libration itself, which is bounded — at least 
while we are far from the libration to rotation transition. Linear stability is not the right tool for this 
mode: we should use orbital stability instead. 

• Eigenvalue 1 diverges for in-plane oscillations for 0 M =54~55 deg., while remaining bounded for higher 
values of 0 M - Resonance rather than de-stabilizing low tension must be at work here. Initial conditions 
are important: initial deviation causes a faster divergence. 

• Eigenvalues 2 and 3 diverge after a number of periods for in-plane oscillations and for the higher 
libration angles; Eigenvalues 4 and 5 remain bounded for the number of oscillations computed (100). 
Again, initial deviation is more critical than initial speed. The fact of being in phase or out of phase with 
the underlying librations seems to matter. 

• Ail eigenmodes except 0 remain bounded for out-of-plane oscillations for the time computed. 

We should look more closely at these points. The behavior of 0 itself, the 0 eigenfunction, will be studied later 
on more in detail. For divergent eigenfunctions 1-3 in-plane, our problem is how long they will take to get out 
of control. We have computed the periods of the oscillations as well, but it is more convenient to display the 
oscillations vs. the true anomaly of the COM orbit, given in number of complete orbits. At the same time, we 
can give a wider variety of initial conditions to further explore the effect of having the oscillations in phase 
and out of phase. This is done by “TOVIC.C” and ‘TOVICNEG.C”. The results show that the magnitude of 
the initial conditions does not matter at all, the equations being linear; the non-linearity and forcing terms, 
that could affect the nature of the solution, are limited to the libration equations and do not carry over to the 
linearized small deviation equations. For all eigenvalues losing stability, the critical case is initial deviation, 
except for eigenvalue 0, which de-stabilizes for initial speed — in phase with the libration. Sign does not affect 
significatively the results. 

The resonance in Eigenvalue 1 for in-plane oscillations is investigated by sweeping more thoroughly the 51- 
61 degree zone. This is done in “TRVS100B.C”, computing 100 oscillations for each eigenvalue and 
independent solution. Results plot oscillation amplitude versus true anomaly of COM orbit, measured in 
number of orbits. The results are shown in Figures 3-7. Figure 3 shows the amplitude of the libration 
computed with the full non-linear equations. El for Eigenvalue 0 in Figure 4 is the variation of the libration 
for initial speed perturbation, while E2 in Figure 5 is the initial deviation variation. 

There is a resonance between 54 and 55 degrees. There is no resonance for out-of-plane oscillations. 
Curiously, there is no loss of stability for higher angles, and for 60 and 61 degrees oscillations remain 
bounded. For Eigenvalues 2 and 3, destabilization follows increasing libration angle, as was noted in 
“TRVSAIOO.C,” “TOVIC.C,” and “TOVICNEG.C.” 

This raises the point of whether there is another resonance at work with these eigenvalues for a higher 
libration angle. The case is not worth pursuing, though, since for these angles the tether goes slack and the 
equations used would no longer be valid. For null tension or even for very low tensions the model used would 
not be valid, and bending stiffness and coiling stresses would have to be considered. Tether points would start 
to move freely until they stray enough to make the tether taut again, causing a rebound. This is a scenario to 
be avoided at all costs. 

Another point to study is the sensitivity of this resonance to variations in the eigenvalue. This will be 
considered after the Ritz-Galerkin method is used. 
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Figure 3 Libration amplitude (full non-linear equations) vs number of orbital periods, for 
different initial libration amplitudes (0 M a X )- 
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Figure 4 In-Plane Initial Speed, 100 Oscillations (El): Amplitude vs number of orbits for 
different initial libration amplitude (0 M ) and for the first six eigenvalues. 
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Figure 5 In-Plane Initial Deviation Oscillations (E2): Amplitud e vs number of orbital 
periods for different initial libration amplitude (6 M ) and for the first six 
eigenvalues. 
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Figure 6 Out-of-Plane Initial Deviation Oscillations (Zl): Amplitude vs number of 
orbital periods for different initial libration amplitude (6 M ) and for the first six 
eigenvalues. 
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Figure 7 Out-of-Plane Initial Deviation Oscillations (Z2): Amplitude vs number of orbital 
periods for different initial libration amplitude (0 M ) and for the first six 
eigenvalues. 
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Ritz-Galerkin Method with Trigonometric Functions 

We will Choose a set of simple orthogonal functions meeting the boundary conditions for the eigenvalue 
spatial problem, and try to minimize the residuals [9,12]: 

■ (36) 
i=l 

with a suitable number of terms k. As a first trial, we will choose trigonometric functions: 

4> i = K, cos(27tiv) + sin(27civ) (37) 

Where the constants K; will be determined so that the boundary conditions (19) are met. Conditions (23) to 
(25) do not make sense here since we are not using separation of variables. These turned out to be the same 
for functions meeting the differential equation. Curiously, for functions (37), both boundary conditions are the 
same as well, but for different reasons: the integral in (19) will be zero for any periodic function integrated 
over a multiple of the period, and for trigonometric functions the second derivative has the same zeros as the 
function itself. Taking four terms for the trial solution, the constants are: 


i 

1 

2 

3 

4 

K 

-0.5955 

-1.8459 

24.5599 

1.5335 


Table 2 

Substituting the trial function in the differential equation, we obtain: 

t 

\ , - H(xK 1 - G(t) j[(b - 1) 2 - V 2 ] ;; } =0 (38) 

” Ha^j --y j[(b-l) 2 - v 2 Ja i <))' , -2va i <t)'}j = R = 0 (39) 

We will use Galerkin weighing to minimize the residual: 

£(}) i R(a j ,(l) J )dv = 0 ) i, j = 1.--4 (40) 


For which we will have to compute the integrals: J^b-l) 2 -v 2 ] <J)" ; which we can put 

in matrix form: 


C 

4>> 



<t>4 

<t>. 

0.6773 

0 

0 

0 

<(>2 

0 

0.3395 

0 

0 

<(>3 

0 

0 

-0.5570 

0 


0 

0 

0 

0.2232 


Table 3 


CP 

vd>,’ 


vfe’ 

V<t> 4 ’ 

<t>. 

-0.3387 

1.6290 

-10.7285 

0.5682 


-1.6290 

-1.1019 

-61.9249 

-2.5623 


10.7285 

61.9249 

-151.0478 

77.1426 

<t>4 

-0.5682 

2.5623 

-77.1426 

-0.83789 


Table 4 
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cs 

f ...Mh” 

l..w 

Ulfc" 

[...]*«” 


-15.8799 

8.6882 

-48.2787 

-2.4244 

_k_ 

2.1720 

-203.3531 

-445.8599 

.13.6656 

4b 

-5.3643 

-198.1599 

-62532.1924 

7053039 


-0.1515 

-3.4164 

396.7334 

-616.0203 


Table 5 


Defining the vector: a = [a, a 2 a 3 a 4 J , we can write the matrix equation: 

C(a-Ha)-— (CS-2CP)a = 0 (41) 

a = [ffl+GC-'(CS-2CP)/2]a = (HI+GA)a (42) 

Where I is the unit matrix and A is: 


-11.2219 

4.0083 

-19.7986 

-0.9507 

1.2320 

-45.6377 

-73.0591 

-1.9378 

-0.0444 

-0.5329 

-102.9973 

9.9119 

-0.3843 

-2.5484 

164.4063 

-186.3008 


Table 6 


Elements in the main diagonal are big and negative, which is good. The trial functions are not eigenvalues, 
but <(>i and O 2 come close to eigenfunctions 3 and 5, as can be seen from the proximity of the diagonal terms. 
Now we have to integrate the time equations for a, but this has to be done numerically: G and H are functions 
of x and 0, which can be obtained as a function of x only through numerical integration. Again, in order to 
avoid integrating Kepler’s equation at every step, we shift to true anomaly as variable. The equations will be: 

„ _ 2esin_v^ + [H( v )i + G(v)Aja (43) 

Where H and G are now those in (31). There is a great similarity with the separation of variables solution. If 
A were just diagonal with the eigenvalues as diagonal terms, the equations would be exactly the same. As in 
the previous case, the integration is made by taking two independent particular solutions and running them 
for a number of orbits. We will take only the two first functions, ai and a 2 , and integrate ala, alb, a2a, a2b 
and 0 for a number of orbits. A is the initial speed solution, and B the initial deviation solution. This is done 
by “RITZTRIG.C”. The results are shown in Figure 8. All solutions show the behavior of the most unstable 
eigenmodes of the previous section. This is to be expected, since trigonometric functions are not 
eigenfunctions of the problem, and each one of them includes a mixture of all modes. The coefficient 
computations can be checked in the MapleV script in Appendix A. 
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Figure 8 In-plane oscillations for Ritz method with trigonometic functions: A1 (first 
mode) and A2 (second mode); oscillation amplitude vs number of orbits. 
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Ritz-Galerkin Method with Legendre Polynomials 

The results of the previous section do not show much detail, since the trial functions are not eigenfunctions 
and therefore each one includes a mixture of all modes. If we want a more detailed output, the trial functions 
should be as close as possible to the real eigenfunctions. This is easily done using the approximate analytical 
solution introduced in (27). Rearranging it and reducing the number of repeated arbitrary constants, we 
obtain: 


A(x) = C 0 


P„ (x) 1 + eK + 2X t eJ o P n (u)Q „ (u)du j + Q „ (x)e 


4(^ 0 -i)r(i+f 


-2^J o X P» : 


du 


(44) 

Where we had made the change v=x(l-b). Since this solution is an approximation, it will fulfill the 
differential equation and boundary conditions (26) only up to 0(e). The value of the constant K cannot be 
determined from the 0(1) and 0(e) problems, we would need to go up to 0(e 2 ), which is already 
unmanageable. But this is not a problem when we use Ritz method: we are not dealing with the Sturm- 
Liouville equations (22), but with the general problem (17) with boundary conditions (19). We can then 
choose K so that the boundary conditions are met. Still, seeing how close our trial solutions come to meeting 
equation (26) tells us how close we are to the real eigenfunctions. Going back to the variable v, our trial 
functions will be: 

Sl( v .'O sss X C i( X ) A i( v ) (45) 

1 = 1 

As in the previous section, only in-plane oscillations will be considered, since out-of-plane oscillations were 
all well behaved — except the fist eigenvalue, which is not a real problem and will be treated separately. 

We will choose IQ so that the conditions (19) are met. This is done with a symbolic manipulator, MapleV. 
The results are collected in the MapleV script in Annex B, showing also how closely the spatial differential 
equation and the different boundary conditions are met. In the computation of Qo(x) we use a recurrence 
relation from [8, § 8.6.19]. The values of K selected to minimize conditions (19) are: Ki=0, K 3 =1.7418, 
K 5 =0.6094. Thus our trial functions are: 


C 1 (t. v) = X c i ( T )<!> i C v) = X c i ( T ) A i f "Tt] • • = 1.3.5 

i=l i=l \1 VJ 


(46) 


Substituting this in (18), we obtain: 
£{c,A 1 -H(T)c l A 1 -M[ (b ..,)-v«]^ r £j 


by 


A „ + G(T)v C | i A , | = r( c . t Aj> v ) , o (47) 


As we did above, we will minimize the Galerkin-weighed residuals: 

b 

jA i R(c j ,A J ,v)dv=0 (48) 

b-1 

For this we will have to compute 27 integrals. It will be simpler to return to the variable x: 

e e e 

(1 _ b)J AjAjdx; (1-b) 2 JxA'Ajdx; (l-b) 3 J(l-x 2 )A;A j dx; i. j = 1,3,5 (49) 

-I -i -1 

This can be put in matrix form defining the vector c and the matrixes A, AP, and AS for the time functions, 
the integrals of function products, the integrals of functions times first derivatives, and the integrals of 
functions times second derivatives, respectively. The integrals are computed with MapleV, as shown in 
Appendix C, and their values are: 
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A 

$1 


... + 


0.3336 

0.01171 

-0.005983 

fa 


0.2226 

-0.0006923 

<f>3 



0.1260 


Table 7 


AP 

x<t>r 

X<k’ 

xfc* 

<f>i 

0.3336 

1.163 

1.097 

02 

0.01171 

0.5938 

1.280 


-0.005983 

0.008283 

0.5253 


Table 8 


AS 

[l-x 2 ]<t>,” 

[1-X 2 ](|)2” 

fl-x’lfc” 

<M 

0 

2.186 

2375 

02 

0 

-1.267 

2.540 

<t>3 

0 

0.04185 

-2381 


Table 9 


A(£-Hc)-yASc + GAPc = 0 


c = 


HT-t-GA' 1 ! — - AP 


c = (HI+GB)c 


Where I is the unit matrix and the value of B, again from MapleV, is: 


(50) 

(51) 


-1 

-0.1663 

-.0305 

0 

-5.5120 

-0.0829 

0 

-0.0693 

-13.6147 


Table 10 


Since the function we used was an approximation, and the boundary conditions were approximately met too, 
there is no point in keeping too many digits in B. We can see that the matrix is decoupled to nearly order 
O(e), and the terms in the main diagonal are very close to the analytical eigenvalues in Table 1 . 

If the boundary condition for the first derivative had been taken into account in the optimization of K„ the 
diagonal terms would be closer to the eigenvalues and the modes would be more neatly separated; but the 
basic assumption of Ritz method — that trial functions meet boundary conditions — would not be as accurately 
met and the results would be less reliable. 

For the numerical integration we turn again to true anomaly instead of non-dimensional time. The equations 
to be run are: 

* \ 3cos0sin0 

0 + 1 

. K (52) 

c= 2eSinV c + (HI + GB)c 

K 

With initial condition (35) for 0, and a pair of independent solutions for c, which is a linear system: initial 
unit speed and zero deviation (a), and initial unit deviation and zero speed (b). This means integrating a total 
of 14 first-order equations. Solutions will be cla, clb, c3a, c3b, c5a, and c5b. The code for this is in 
“RITZLEG.C ” cl corresponds to the eigenvalue we earlier called 0, which is the rectilinear shape; c3 
corresponds to eigenvalue 1 of previous figures, and c5 to eigenvalue 2. 

The results are shown in Figure 9. All solutions grow exponentially after a number of orbits. This is similar to 
the previous case, in which trigonometric functions encompass all eigenmodes and diverge quickly, while for 
the approximate eigenfunctions there were clearly stable and unstable ones. In this case the onset of instability 


0 = 


2esin v 
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can be traced to the coupling terms in matrix B. Although small — say of order e — their effect would be felt 
after a time of order 1/e. Thus cl will feel the influence of c3 after 1/0.16, or about two orbits. This is not 
noticeable because the second eigenvalue takes about 30 orbits to diverge. It would feel the influence of c5 
after 1/0.03, or about 10 orbits, and so it diverges at that time. The initial deviation solution diverges earlier 
because the corresponding c5b does so as well. Finally, c3 should diverge after 1/0.08, or 6 orbits, and indeed 
it does so. 

Thus the off-diagonal terms of matrix B couple the modes and mask each one’s proper behavior. Short of an 
exact eigenmode, we can simply omit these terms and keep the main diagonal, or the approximate 
eigenmodes. This is done in “RITZLEG2.C.” The results are shown in Figure 10. Each mode recovers its own 
personality, and we notice even the decay of the initial deviation clb solution (eigenvalue 0 previously) for 42 
degrees, that was present in the “TRVSA100.C” results (not shown in Fig. 5 because the “TRVSB100.C” 
results scan only 51-61 degrees). Of the 0 M values computed, this is the closest to 45 degrees, where the 
maximum of the restoring term 3sin0cos0 is experienced. 
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Figure 9 In-plane oscillations for Ritz method with Legendre functions: 3 first modes; 
oscillation amplitude vs number of orbits. 
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Numerical Determination of Eigenvalues 

The approximate solution (27) to the Sturm-Liouville problem (22) cannot go beyond the 0(e) expansion, 
since the equations become too complex. There is always the possibility of numerically solving the differential 
equation. A tabulated solution is not easy to handle, but what we need for the time equations is only the 
eigenvalue. We set then to solve equation (22) with boundary conditions (23), (24), or (25) which, as we saw 
before, are exactly the same. (23) or (25) are easy to impose, while (24) can only be used after the integration 
to check the accuracy of the solution. 

The only boundary condition at the free end is that the tension be zero, which causes the coefficient of A” to 
be zero. A”(b-1) is thus an indetermination of type 0/0, and cannot be used to start a run. What we know is 
that it cannot be infinite, since it is related to the curvature of the tether. Therefore, the rest of the differential 
equation must be zero as well, which provides a boundary condition to check proximity to the real eigenvalue. 
We can thus start integrating at v=b with boundary condition (22) and an arbitrary value for A(b). Since 
eigenfunctions are defined except for a multiplicative constant, we may normalize the value of A so that, for 
the linear solution, A(b-l)=l, and thus A(b)=b/(l-b). 

At the free end we will require that A” be finite, that is, 

-(1 - b) A'(l - b) ■ + AA(1 - b) = 0 (53) 

which is the same condition imposed on the end mass, (25). The functions to integrate are: 


i 

yBl 

dy[i] 

Init. Cond. 

1 

1 

y[2] 

0 

2 

A 

y[3] 

b/( 1 -b) 

3 

A’ 

(2vy[3] - 2Xy[2])/[(b - 1) 2 - v 2 ] 

A/(l-b) 

~T 

X 

0 

k(k+l)/2 (1,3,5...) 


Table 11 


With this we start a shooting method integrating routine in order to determine X. The residual of (53) is taken 
as a function of X whose roots we try to find with a combined shooting method/Newton-Raphson algorithm 
[6, pp. 756ff]. This is done calculating the Jacobian (which in this case is just the derivative, since we have 
only one variable): 

- *fc>l _ X , .!,_«£■> (54, 

+AA. -» R^+AX)] A X J 

We use as seed value the approximate eigenvalues found by asymptotic expansion, (27). y[l], which is 
condition (24), can be used as a further check of our approximation. 

Still another simplification can be made by a change in variables: x=-v, which avoids integrating backwards 
and checking all integration routines for the sign of h. The form of the differential equation does not change 
at all: 

[(b-1) 2 -x 2 ] A"-2xA' + 2^A = 0 (55) 

And the boundary conditions remain unchanged, except that they are imposed at -b; the integrating interval 
becomes [-b:l-b]. The functions to be integrated now are: 


i 

y[il 

dy[i] 

Init. Cond. 

1 

j 

y[2] 

0 


A 

y[3] 

-b/( 1 -b) 


A’ 

(2xy[3]-2Xy[2])/[(b-l) 2 -x 2 ] 

X/(l-b) 

4 

X 

0 

k(k+l)/2 (1,3,5...) 


Table 12 
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There is an additional problem: the equation is singular in x=l-b, the free end. We cannot integrate all the 
way to the end, because A” will be unbounded in all iterations short of the exact eigenvalue. To sidestep this 
singularity, we will stop within a distance 5 of the end. 5 is chosen as ~10' 8 so that it is much smaller than the 
approximation we seek (about e 3 , with e~0.085), while still big enough that division by 8 2 will not cause 
overflow for double precision arithmetic. Otherwise, the error in each integration step goes out of hand as we 
approach the end. 

This is done by “EIGVAL.C,” using routines taken from [6] like fdjac.c, fmin.c, lnsrch.c, while 
others have been modified to fit this problem, like newt.c, odeint.c, and shoot. c. The first six 
eigenvalues are computed with this program. 


X 

Numerical 

Analytical 

Ritz-Legendre 1 

Ritz-Legendre 2 

xo 

1.000000 

1.000 

1.000 

1.000 

A.1 

5.439627 

5.386 

5.512 

5.477 

X2 

13.453447 

13.312 

13.614 

13.558 

X3 

25.031209 

24.764 

- 


X4 

40.172217 

39.740 

- 

_ 

X5 

58.875815 

58.240 

- 

- 


Table 13 


The approximate analytical values come short of the numerical values, while the diagonal terms of matrix B 
in the Ritz method with Legendre polynomials outstrip them. The difference is greater when only the Ritz 
method-specific boundary conditions (19) are used to optimize the free parameter K (column 1). If condition 
(25) is included in the optimization (column 2), the diagonal values are closer to the numerical eigenvalues, 
but the Ritz method results are not as reliable. 

These numerical eigenvalues are used to compute again the transversal oscillations. This is done by 
‘TRVC100.C.” 

The reason for repeating the same computations every time we find some small change in the eigenvalues is 
that the stability of the modes may depend sensitively on the eigenvalues. For a circular orbit, the equations of 
the transversal oscillations become Mathieu equations [3, p. 386], whose stability areas are well known. The 
unstable areas increase for larger libration angles, so that a change in the eigenvalue may, depending on 
where we stay, push a mode in or out of the unstable area for a certain libration angle, or could raise or lower 
the libration angle at which a mode loses stability. 

When we are in an elliptical orbit, more terms enter the equations and we have the additional resonance with 
the orbital frequency; but for small eccentricities we expect some similarity with the Mathieu pattern. A new 
run is thus well justified. The results are collected in Figures 1 1-14. A gain, remarkable points are: 

• The initial speed El solution shows that de-stabilization of eigenvalue 1 is pushed up to near 55 degrees., 
while in Fig. 4 it was closer to 54 degrees. Other eigenvalues keep the same behavior, with small 
variations in amplitude indicating that we are getting closer or farther from the unstable areas. 

• The initial deviation E2 shows the same shift towards 55 degree de-stabilization as El. Eigenvalue 5 
takes longer to diverge. 

• Initial speed out-of-plane Z1 shows no change with respect to Fig. 6, apart from eigenvalue 3 that shows 
a beginning of de-stabilization which was not present in any of the out-of-plane eigenvalues except 0. 

• Initial deviation Z2 shows two unstability zones for eigenvalue 0, for 60 degrees and larger, and for 53 
degrees; the latter just started to show in Zl. It seems that, as 0 M grows, we enter an unstable Mathieu 
zone, go out of it as it twists left, and finally enter the second unstable zone for 60 degrees and larger. 
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Figure 14 Numerical eigenvalues, out-of-plane initial deviation (Z2), 100 Oscillations: 
Amplitude vs number of orbits for different 0 M and for the first six eigenvalues. 
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Unstable Modes 

The previous analysis of transverse oscillations has singled out some unstable modes and de-stabilizing 
factors that need further study. At the same time we should check the limits of our assumptions, like the 
decoupling of transverse oscillations. Points to be studied are: 

• Relation of tension and stability 

• Zero mode divergence 

• Out-of-plane oscillation coupling 

• Libration to rotation transition 


Tension and Stability 

For transverse oscillations, tension is the main stabilization factor, providing the restoring force for the 
oscillations. As tension becomes smaller, oscillations tend to increase in amplitude. Pendular libration and 
tension are governed by equations (6): 


0 - 


2e sin v 


sin 0 cos 0 = 0 


F = 


or 

2n* 


(e+i)+- 

\ / K 

(0 + 1) 2 - - (l - 3 cos 2 e)l[(b - 1) 2 - V 2 ] = G( v)P( V) 



Non-dimensional Tether Length 


P(v) reaches its maximum at the COM, and is zero at the free 
end. The tension distribution keeps its shape at all times, its 
magnitude governed by G(v). When G becomes zero, the tether 
goes slack at all points at the same time, and from that moment 
on the equations used are no longer valid. The tension time 
dependence, G, appears in the restoring term of transverse 
oscillation equations (28); 


Z + 


r co 2 ^ 

XG(x) + -^ t 

ten j 


Z = 0 


E + (AG('t) - H(t))E = 0 


Figure 15 Tension distribution It is thus interesting to see when the tension and the restoring 

terms go to zero. This is done in “THETA20.C,” computing 0, 
G, and H for a number of orbits, the results are shown in Figures 16 and 17. 

Theta remains bounded, and the tension is positive up to 0 M =63 degrees. G goes to zero between 63 and 64 
degrees. There is a zone where H is bigger than G, so that the restoring force is negative for low eigenvalues, 
specifically for X=1 (eigenvalue 0). Since X increases fast, for eigenvalues 1 and above, the value of H is not a 
problem except when G is very close to zero, that is, for large libration amplitudes. For angles above 63 
degrees, G has negative zones, so that larger values of X mean broader negative restoring force zones. 
Out-of-plane oscillations are more stable because of the positive go 2 term. This is just part of the picture. 
Resonance with libration and orbital frequency has to be taken into account as well. We will look into these in 
the following. 
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Figure 16Libration angle for different initial libration amplitudes. 
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Instability of Eigenvalue 0 


Eigenvalue 0 is unstable for initial speed in-plane and for all out-of-plane transverse oscillations. Initial 
deviation in-plane oscillations remain bounded (at least for the number of orbits computed), and even shows a 
decrease for 42 degrees (cf. Fig. 10). 

Eigenvalue 0 is a rectilinear shape, and corresponds to a variation in the initial conditions of the pendular 
motion. From the equations of pendular motion (3), introducing non-dimensional time x, we obtain: 


0 - 


<P + 


2esinvco 2 


ten 


0 + — 


- 2cptan <p 


)j + “' 


3co 


Kn 


cos0sin0 = O 


3oo 2 


Kn 


cos 2 0 


(56) 


sincpcoscp = 0 


For libration in the orbital plane, cp is 0, so that we can consider only 0. Giving an infinitesimal variation to 0 
we get: 


0 + 80 — 


2e sin vco 3co 


cos(0 + 80) sin(0 + 80) = 0 


Kn Kn 

Considering (56) and dropping higher order terms, we get: 

3cd 2 

80 h r [cos 2 0-sin 2 0j80 = 0 

Which is precisely equation (28) for in-plane transverse oscillations E and ^=1. In the 
substitute an infinitesimal deviation 8<p for the <p=0 out-of-plane libration solution, we get: 


8cp + 




+ — -cos 2 0 |5<p = 0 
Kn 2 


(57) 

(59) 

same way, if we 


(60) 


Which is precisely equation (28) for Z and h=\. Therefore the increase in Eigenvalue 0 amplitudes is not 
worrisome. There is a monotonous linear increase in El with subharmonic modulation, as shown in Figures 
4, 10 and 11. Z1 and Z2 show an exponential growth. This corresponds to orbital stability. The evolutions of 
0 and cp are non-linear oscillations, with the period depending on the amplitude (plus the forcing term due to 
the elliptical orbit, which is small). A change in initial conditions means a change in period, so that the small 
initial deviation will increase continuously, but not go beyond the amplitude of the new oscillation, which 
remains bounded as long as we are far from the potential crest. We will consider later on at what initial 
conditions inertial orbital period forcing moves the libration into rotation. 


Z1 and Z2 are deviations from a zero amplitude non-linear 
oscillation, in which case the linearization is quite accurate; 
however, the restoring term varies with time, which may cause 
resonance. This accounts for the exponential amplitude growth 
of Z1 and Z2. 

El and E2 behave differently, the latter remaining bounded for 
the time interval computed (50-100 orbits). The libration 
movement is essentially a gravity gradient pendulum (soft 
spring) with a small forcing term due to ellipticity of the orbit. 
The small forcing term will have an accumulated effect when 
Figure 18 Gravity gradient pendulum there is resonance, but for a short span can be left out. This 
potential leaves a conservative force field with a potential function, as 

shown in Figure 18. An infinitesimal deviation when the 
pendulum crosses the vertical means moving along a nearly horizontal potential curve (deviation line in Fig. 
18). This will add to the total energy an amount one order of magnitude smaller. For libration amplitudes of 
about 60 degrees a speed differential means moving up a potential curve with a high slope (speed line in Fig. 
18). This will add to the total energy something of its own order, thus affecting the amplitude and the period. 
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Figure 19 Eigenvalue 0 and the difference between two librations. 

This can be checked by comparing the Eigenvalue 0 deviation for 60 degrees libration with the difference 
between the basic libration (60 degrees) and the new libration obtained by adding the initial conditions of the 
Eigenvalue 0. We will use for initial speed the speed difference at zero crossing between a 60 degree libration 
and a 61 degree libration, while for initial deviation we will just push the 60 degree libration one degree 
away. Figure 19 shows the agreement between the libration increment and the differential. 

The initial conditions for initial speed and deviation are: 

e,(0) = 0 0 1 (O) = V3sin6O 0,(0) = 0 0,(0) = V3 sin60 

0 2 (O) = 0 0 2 (O) = V3sin61 0 2 (O) = 1 0 2 (O) = V3sin60 

80(0) = 0 50(0) = 0 2 (O)- 0,(0) 50(0) = 1 50(0) = 0 

As shown in Fig. 19, the agreement is quite good, given the fact that eccentricity has been assumed to remain 
constant, while in fact it depends on the initial libration conditions. The difference, however, will be of a 
smaller order. Thus the divergence in Eigenvalue 0 is not a problem. 


Out-of-Plane Oscillations and Coupling 

We will now consider the out-of-plane libration. Z1 and Z2 lose stability faster than El for Eigenvalue 0, 
which is the linearized libration. We will compare these solutions with the full non-linear libration for small 
amplitudes. Changing variables in equations (56) and (60) from x to true anomaly v we obtain: 

q = ^ es * n v (6 + l) + 2<ptan cp^0 + 1) - — cos 0 sin 0 

K K 

2esinv . 

9 = (p- 

K 


(0 + l) + — cos 2 9 


sin<pcos<p 


(61) 


.... 2esin v c . 
ocp = ocp — 




cos 2 e 


5<p 


Still, in order to make the initial conditions comparable, we have to relate the initial speed and deviation. For 
this we assume a small oscillation whose simplified equation would be: (p + co 2 <p = 0 . An initial deviation 
solution of amplitude A would have a speed Aco when crossing the orbital plane. We compute co from the 
value of the &p coefficient for 0=0 and 0 M =6O degrees: 


co 




(62) 


We will consider that <p is small, so that the tp term in 0 equation will be small as well. For the coupling term 

I <p = 

to be smaller than the velocity term (order e), tp should be: (p tan <p = e => < r . 

[tantp = cp = ve 
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Figure 21 Influence of <p on the orbital plane libration for 0 M =6O deg. 


We will integrate the full libration equations for 60 degree in-plane libration, and out of plane oscillations 
ranging from 0 to 10 degree amplitude. The initial deviation will correspond to full amplitude, and initial 
speed to the plane crossing speed for that amplitude according to (62). 

The results are shown in Figure 21. In both cases, there is an effect over the amplitude and the period of the 
libration, which becomes smaller as tp increases. For the time — 20 orbits — and the out-of-plane oscillations 
computed — 0 to 9 degrees — the effect is not enough to push the libration into rotation. It is enough, however, 
to disrupt significatively the libration, so that the computations made to maximize the momentum transfer to 
the satellite would no longer be valid. 

Consequently, out-of-plane oscillations should not be allowed to increase beyond a few degrees of amplitude. 
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Libration to Rotation Transition 

In what we have been studying, the assumption was made that the tether librates without going into a 
rotation. However, the gravity gradient pendulum is a sort of soft spring, so that at large amplitudes we must 
be wary of going over the potential crest. The inertial term acts as a periodic forcing, which may force the 
libration into rotation when we are close to the limit, or have a cumulative effect through resonance. We will 
see what are the libration angles and speeds leading to rotation. Even though it may happen at angles that 
will never be reached — tension becomes zero between 63 and 64 degrees — it is good to know how safe we are 
at 60 degrees. 
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Figure 22 Libration phase maps for initial deviation. 20 orbits. 

Resonance is an important factor. Initial deviation means cutting the tether at extreme amplitude and no 
speed. This sends the TS into a low orbit with a smaller eccentricity. Still, there is a transition to rotation for 
66 and 67 degrees, while for lower and higher angles there is no transition, as shown in Figure 22. This is an 
effect of resonance with the forcing term, which increases the energy of the libration. The initial deviation 
condition is far from our case, since we will be cutting the tether at vertical crossing, which is a pure initial 
speed condition. 

The case we are most interested in is the initial speed libration. Here, the transition is more regular. As shown 
in Figure 23, amplitude keeps growing, but there is no transition to rotation within the interesting libration 
angles. 
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Figure 23Libration phase maps for initial speed, 50 orbits. 
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DECAY MINIMIZATION 

The purpose of the swinging tether, double cut procedure is to maximize the useful life of the satellite. We 
seek the combination of parameters that will yield the maximum orbital lifetime. For this, we will seek to 
minimize the initial orbital decay rate as a function of the deployment free parameters. 

Assuming that the Shuttle is in a circular orbit, the deployment strategy leaves three free parameters: 

• Libration amplitude, which depends on the tether deployment strategy. 

• First cut timing, or gravity gradient pendulum libration angle at cut time. 

• Second cut timing, or libration angle of the TS in the elliptical orbit frame. 

The Shuttle will be in a circular orbit at an altitude of 297 Km, with 57 degrees inclination. The tether is 20 
Km long, with a linear density of 3.3 1 O' 4 Kg/m, and the end mass is m A =32 Kg. We assume the earth to be 
spherical, its radius equal to the equatorial radius. This is a conservative assumption, since the orbital altitude 
will be greater than the equatorial altitude for most of the period, with a correspondingly lower air density 
and lower decay rate. Thus the circular orbit radius will be Ro=6675.16 Km, and its angular rate will be coo- 
After deployment, the TS librates as a gravity gradient pendulum governed by equation (5). Since e is zero, 
the equation admits an energy integral: 

9 2 + 3 a>l sin 2 0 = 3©o sin 2 0 M (63) 

Where 0 M is the maximum amplitude of the gravity gradient pendulum libration. We keep the same reference 
frame we have adopted previously, with Ox as the local vertical and Oz in the orbital speed direction. 8 will 
be the libration angle at cut time. Lf will be the final deployed length of the tether; tf will be the distance 
from the TS center of mass to the free end, that is: ^=oLf, with a=l-b as defined in (4). 

We will compute the parameters of the different stages of the TS motion as a function of the deployment free 
parameters. 


Elliptical Orbit of the COM 

When the tether is cut at the orbiter end, its center of mass, which we will call G, enters an elliptical orbit [2] 
with the initial conditions: 


(64) 


The initial libration angular speed can be obtained from the energy equation (63): 0 O = <n 0 u , with: 
u = >/3.y/sin 2 0 M - sin 2 8 (65) 

which is different from the parameter u used in [2], Since the tether is much shorter than the orbital radius, 

f 

we can use an expansion in powers of a small parameter — = e « 1 , leading to: 

R„ 


Ro =• 

Rq 4* COS 5 

0 

II 

O o 

> 

(o) 0 + 0 O )^j? sin 8 
0 


sin 5 


^0^0 4" (CO 0 4" 0q)^f ^ 


Ro 


1 + ecosS 


-(1 + u)esin8 

0 

Y © 

04 

o 

3 

II 

O o 

> 

0 

esinS 


1 + (1 + u)ecos8 


The orbital parameters can thus be expressed as follows: 

r 7-r- — — 7 (cos 2 8 + lY9 + 18u + 12u 2 +2u 3 ) , 

e = v3cos 2 8 (u + l)(u + 3) + u e + cos 8 - — A 'c 2 i 

2^/3 cos 2 8 (u + l)(u + 3) + u 2 

a = R 0 jl + 2 cos S(u + 2 )e + [cos 2 8(l3 + 16u + 4u 2 )+ 2 + 2 u + u 2 j e 2 +. ..J 
n = co 0 jl -3cosS(2 + u)e--^J(u 2 + 2 u + 2 )-cos 2 8 (u 2 + 4 u + 7)j e 2 +...| 


( 66 ) 

(67) 

( 68 ) 
(69) 
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( 70 ) 


Ah p = a(l-e)-R 0 = R 0 [“V^ cos2 S(u + l)(u + 3) + u 2 + 2 cos 8(u + 2)je+. . . 
Ah a = a(l + e)-R 0 = R 0 |^/3cos 2 5 (u + 1 )(u + 3) + u 2 +2cx>s5(u + 2)je+... 


For small 5, these expressions reduce to the equations derived in [1, p. 4-5] and [2, pp. 5-6]. Since we are 
precisely trying to study the influence of 5 in decay rate, we cannot assume at the beginning that it is small. 

r a e 

We still need the argument of the perigee, which can be obtained from: sin v. = — — • j , which gives; 

M-H 


u sin 6 ^ e cos 2 S(2u" + 6u + 3) - (4u + 3) 

v = - arcsin| -usino r — £+... 

2v 2 


v 


( 71 ) 


where: v = ^3 cos 2 8(u+3)(u + l) + u (72) 

With the expressions above, the orbit of the COM is perfectly determined. The details of the mathematical 
derivations and simplifications are in Appendix E. 


Libration of the Tethered System in an Elliptical Orbit 



Figure 24 Reference frames 
for the first cut. 


After the first cut, the TS continues its libration given by eqs. 
(5) while its COM follows an elliptical orbit. A new reference 
frame is erected, Gxyz, with G at the TS COM, Gx along the 
local vertical towards zenith, Gz in the orbital plane and Gy 
normal to the other two. This is the same reference frame 
defined in page 3. As we saw before, equations (5) — or 
equations (6) when we take true anomaly as variable — have to 
be integrated numerically. For this, the initial conditions must 
be determined. 

Referring the index 1 to an earth-centered fixed frame, 2 to 
the TS, 0 to the COM orbital frame, and 3 to the circular orbit 
frame, we seek 0^0=023+0)3 i-g>oi (cf. page 9 ). 

Besides, we will temporarily distinguish between: 

• 0= angle between the rigid tether and the local vertical at 
the Orbiter. 

• \\f= angle between the rigid tether and the local vertical at 
the COM. 

Later on both will be called 0, but initially we have to 
distinguish between them. At cut time we have therefore: 


Vo = 00 ■ 


c* sin 5 r. £ sin § ^ • c 2 o . ^ 

. v n = 8 - arctan — — 0 - arctan r~o-ssino + £ cososino (73) 

R 0 +^ F cos 8 1 + ecosS 


For the initial angular speed, we use (33)-(35), taking into account that the parameters u and £ from [2] used 
there are different from what we are using here. With the new u and £, and substituting (67) and (69), we 
obtain: 


v = n { ' +CC ^) =0) 0 [i+[ Co + C| cosv + 


c, cos V 


]e+...} 


(74) 


with: 

c () = -3 cos 8(u + 2) + e[6 cos 2 8(2 + u) 2 - 3(u + 1)] 

< c, = 2v-ecos5[(2u + l)(2u 2 -9) + (16u 3 +96u 2 + 180u+99)cos 2 8]/ v (75) 

c, = ev 2 
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In order to obtain the initial angular rate CQoi, we only have to substitute the true anomaly at cut time from 
( 71 ) in the angular rate expression ( 74 ): 

. „ (l + ecosV ;) 2 r , , , 

= v(v i ) = co— -375— = co 0 [l + ucos5e-u(2cos 2 5-l)e 2 +...l (76) 

(!-e) 


tn„ 


From (63) we obtain. co 22 — 0 O — to 0 u ; (Oj, — co 0 — ■ Adding all the terms, we obtain the initial 

conditions: 


("cit") =< ° 20 = CD o u [ 1-cos 5 e + ( 2cos2 8-l)e 2 +...] 

=-^ 2L = u|l-cos8(l + u)e + [cos 2 8(u + 2)-l](u + l)£ 2 +...j 


(77) 


We can now integrate the libration equation. We return to using 0 for the libration angle, once the circular 
orbit frame has disappeared and there is no possibility of confusion. As usual, we adopt v as the integration 
variable. In order to obtain the TS angular speed needed in the momentum transfer computations, we will 
have to multiply times the orbital angular speed. We integrate numerically equation (6) with initial conditions 
(73) and (77), up to the true anomaly at the second cut. 


Final Orbit of the Satellite 

The computation process we are following is: 

• Circular Shuttle orbit: Ro, cob- 

• Elliptical orbit of the COM after the first cut: a(8,0 M ), e(8,0 M ), Vj(S,0 M ), co(8,0m), 0(v), 0’(v). 

• Elliptical orbit of the satellite after the second cut: a s (S,0 M ,v c ), e s (8,0 M ,v c ), Vk(8,0 m ,v c ). 

Where v c is the true anomaly of the COM orbit at the time of the second cut. In order to determine the final 
orbit parameters, we will obtain first the position and speed of the end mass (A) as a function of v and the 
COM orbital parameters. 

r A =r°+GA = [R + GAcos0,O,GAsin0j (78) 


Were: R = 


(1-a) 


, and GA = L F - E F = 

1 + ecosic a 


^ F = peR 0 ; with P=0.0934937, If we compare terms: 


8=0.0028 (5=0.0934 

e 2 =0. 0000822 p 2 =0.00874~e 

Consequently, we only need to keep terms in £, p, p 2 , and eP, if we are disregarding terms of order 0(£ 2 ). 
Actually, there should not be any need for analytical computations at this stage: once the libration has to be 
computed numerically, we can simply compute speeds and positions at each requested point and apply 
standard formulas to obtain the orbital parameters. However, the problem is stiff due to the magnitude 
difference between a, Ro, and a-Ro. We will perfect the algebra in order to compute only the increments, 
which yield reasonable results without resorting to double precision. From the previous section we have: 

a = R 0 (l + a 1 e) co = (o 0 (l +(*>,£) e = e,£ v = <o 0 Jl + (c 0 + c, cosv + c 2 cos 2 v)J 

where the perturbations are known and computed. If we desire greater precision, the £ 2 terms can be included. 
The absolute speed of the end mass is: 


V?, =V A +v A =V A +V° +( 


• co 01 A GA = (0 + l)vp£R 0 
Developing in powers of the small parameters: 


-sin0 


esinv 

0 

aco 

+ T j 

0 

cos© 

Vl-e 2 

1 + ecos v 


(79) 
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l + e(3cos0 + a, -e, cosv) + 0(e 2 ) 


1 + er.x 

o 

04 

1! 

< 

!■« 

0 

- = R o- 

0 


Pe sin 0 


. Er iz . 


v 2[ - R 0 °V 


eje, sin v - (3(0 + 1) sin 0]+-*- 



0 

* — Rq^V 

0 

1 n- e[ej cos v + a ] + C 0 [ +p('d + l)cos0]+--- 


. 1 + £v u. 


( 80 ) 


The details of the simplification can be checked in Appendix F. Keeping only the perturbations of r A and v A t 
we compute the parameters of the final orbit: 
a s = R 0 (l + a sl e) 

with: 


a si =2(v, I +r u ) + e(5vf z +vf x + r, 2 z + 2r 2 + 8r lx v lz ) 


«sl = A ~ 


effjz ~ IQfrx V L ~ 2 v iz r iz + 2v lz r, z v lz - 4v lz r 2 + 2v u r„r u + r lx r, 2 z - 2v u v? x - 2r„ v? x ] 

2A 


(81) 


A = V4v?, + 4r lx v lz + r lx + r 2 + 2r, z v lx + vf x 


Again, the details of operations and simplifications are in Appendix F. The expressions above provide the 
parameters of the final orbit of the satellite, as a function of 5, 9 m, and v c . We can then compute the orbital 
decay rate. 


Orbital Decay Rate 


Speed Determination 

In order to find the height loss due to air drag, we will assume: 

• Spherical earth and atmosphere. 

• Atmosphere rotating at the earth angular speed 

• Exponential density model. 

The diurnal bulge will not be considered at this stage. Night/day changes will not be considered either, nor 
solar activity changes: we are seeking a simple analytical model for a first estimate. 

We introduce an earth centered “fixed” frame (system 1); the orbital frame (system 0), which is actually the 
frame of the polar coordinates in the orbital plane; and the satellite itself (system 2). 

Air drag can be estimated [1,7, 11] as: 


D = -|pv; 0 AC Dl ^ 


(82) 


Where p is the air density, A the average cross section of the satellite, and C D the drag coefficient [1, p. 6-8]. 
In order to compute the relative speed we need the orbital speed of the satellite and the rotational speed of the 
atmosphere, projected onto the orbital frame, u^Uv,!^: 


v;, = 


aco 




t [e sin v,l + e cos v,0 J = [e sin v,l + e cos v,0j 


v s =Q 

M)1 ' 


. r 5 with: 


r = ru = - 


a(l-e 2 ) i 
1 + ecos v 


Q, 


= Q E [sin i sin(co p + v), sin i cos(co p + v), cos i J 
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acoe sin v 


Vl-e 2 

aco(l + e cos v) Q E a(l - e 2 )cos i 


vr 


. e 2 1 + e cos v 

Q E a(l - e 2 )sin i cos(co + v) 


( 83 ) 


1 + ecos v 

where C0p is the argument of perigee from the line of nodes, and Q E the earth’s rotational speed. The drag 
expression includes the square of the speed. The air speed is small with regard to the satellite orbital speed, 
and eccentricity is small as well. We will try a Taylor expansion, but we must check the orders of magnitude 
of the quantities involved: 

• Q e = 7.29210' 5 rad/sec [7, p. 350]. 

• £2 e /cd -0.05 

• e~0.02 

Consequently, we can drop terms of order e 2 ,£2 E 2 , and eQ E The squared speed will be: 


v™ = 


(1 + 2e cos v) 


1 - 2 — cos i(l - e cos v) 
Q 


(84) 


where we have used the approximate orbital rate, Q = — = (1 + 2e cos v) « n(l + 2e cos v) . The last term 

r h 

in (84) can also be dropped, since it is of order eQ E This leaves the same expression as in [7]. 


Energy Loss 


dE 


The variation of mechanical energy can be obtained from the drag power: — = D v , with D and v as 

dt 

y s * - 

shown above. The dot product can be approximated by: j-^j- • v s 21 = [vy , since the difference is of order e 2 . 


This leaves just the product of scalar values: 
dE 


dt 


■ = D v* 


(85) 


Orbital Decay Rate 

Semiaxis and eccentricity variations can be obtained from the equations for variations of elements [7, p. 210]. 
The semiaxis variation can be obtained directly from the energy equation: 

_ m|i da 2a 2 dE 2a 2 _ 

2a dt mp dt mp. 

Substituting the value of v, expanding, and dropping terms of order e 2 , we obtain: 

— = (l + 3ecosv) = — (1 + ecos v) (86) 

dt m D. mn 

The second form is obviously easier for integration. 


Semiaxis Decay in One Orbit 

It is possible to integrate (86) for one orbit, and obtain the decrement of the semiaxis. By minimizing this 
expression, we will maximize the life expectancy of the satellite. For the analytical integration, it is simpler to 
shift to eccentric anomaly, E: 

. r T , r 2 * da dt 

Aa= da= dE (87) 

Jo Jo dt dE v ; 
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dt = 


1-ecosE 


We need to relate v and t with E: 
nt = E-esinE 
ndt = dE(l - e cos E) 
cos E - e 

cos V = , 

1 - e cos E 1 - e cos E 

The air drag must be expressed in terms of E as well: 


1 + e cos v = • 


dE => da = 
1 


2D(l+ecos v) 


da = 


mn‘ 

2D 

rmrdE 


(l-ecosE)dE 


D = - ~ P AC D * 


2 ' . cosE-e 

1 + 2e 

1 - e cos E 


1 - 2 (1 - 3e cos v) cos i 

n 


]■ 


= -^pAC D -j~(l + 2ecosE) 1-2-^-cosi 


(88) 


For the air density, we will use the well known exponential model: p =p 0 e H , where H is the scale height 
which varies with the reference height ho. We only have to give h as a function of E and orbital parameters: 
h n - h _ h 0 -[a(l-ecosE)-R E ] _ h 0 +R E -a + ae cosE 

H H ~ H 

Substituting these expressions into the integral, we obtain: 


Aa = -«. e ^ 

mn 2 h 2 




A (2 it < 

-cosi IJ e H 


(1 + 2e cosE)dE 


(89) 


Where we have considered that Aa and Ae will be much smaller than a and e, which , in turn, can be assumed 
constant for one orbit. Thus the integral admits an analytic solution in terms of Bessel functions [7, p. 212]: 

\f ^0 + Re -a '\ / 


^ a 2rt = "I 271 


ACpPo \ 


( P 2 ' 


m y 

The values of the constants used are: 


„ 2 . 2 

V n h J 


t _ ftc 
1-2 — —cos i 
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( ae ^ _ _ ( ae 

•IhMh 


(90) 


AC d 1 1 m 2 1 iri -« Km 2 

= — = = 10 , from [1, p. 8]. 

3 118 Kg 118 


m 


Kg 


Kg 


• p 0 = p(300Km) = 3.3 ■ 10' 14 -ly = 3.3 • 10~ 2 -pj- , from [12]. 

cm Km 

• H(300Km) = 55.77Km , from [12]. 

Consequently, the value of the first parenthesis is: K=0.1757 10' 8 Km' 2 . 

..2 


■ = a“ 


hp^^E ~ a 
* H 


1 - 2ft cos i = 1 - 0.125 • lO-'V? 


where we have used n E =7.29 10' 5 rad/sec and p=398618Km 3 /sec'\ from [5, pp. 350-51]; and i=57 deg from 
the latest unpublished flight data. 

If we finally introduce the orbit parameters in the main term plus the small-perturbation form obtained in 
(81), we have: 

Aa = 0.1 7(0.4 1 5 + 0.002 la„ + 0.0000023a 2 ) e 0 05 -° 306a - [l 0 (z) + 0.00541, (z) ] (91) 

Where: z = 0.45 • 1 0' 2 (66.75 + 0.1 8a sl ) e sl . 

Again, the computational details are shown in Appendix H. 
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Computation Results 

We expect that the maximum perigee height, which is the main concern for orbital decay, will be attained if 
the tether is cut at apogee, with zero libration and positive libration speed, so that additional momentum is 
transferred to the satellite. We thus compute the orbital decay rate as a function of three parameters: libration 
amplitude before the first cut, 0 M ; libration angle at the first cut, 5; and COM orbit true anomaly at the second 
cut, v c (the second cut time would be an equivalent variable, as they are related by the Kepler’s equation). 
Integrating for a range of values, we obtain: 

• Minimum decay rate always for 8=0, but with small sensitivity for values 8 of up to 9 M /2. 

• Minimum decay rate when the tether is cut at apogee, with the second apogee passage yielding a lower 
rate than the first; small sensitivity to true anomaly as well. 

• Minimum decay rate for 0 M ~6O degrees. 

The lowest value in the range explored is for: 


Aa: 

0.014541 

Km/orbit 

9m: 

60.000000 

deg 

8: 

0.000000 

deg 

v c : 

559.529114 

deg 

0(v c ): 

-6.386499 

deg 

0’(v c ): 

82.325592 

deg/rad v 


The results can be seen in Figs. 25-31. For low libration angles, the cut at the first apogee passage produces a 
minimum decay rate, as seen in Fig. 25. As the libration angle increases, the boost from the tether increases 
as the libration speed tunes its phase with the orbital anomaly, so that the absolute minimum is reached for 60 
degrees (Fig. 27). The libration speed is then positive and maximum (Fig. 28). The apogee and perigee 
heights obtained are shown in Figs. 29-31 for different 0 M . 
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Orbital Decay Rate for Theta Max = 58 deg 
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Orbital Decay Rate for Theta Max = 60 deg 



Figure 27 Decay rate (Km/orbit) for libration angles about the minimum (60 deg). 
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Apogee Height for Theta Max = 58 deg 
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Perigee Height for Theta Max = 59 deg 
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Figure 29 Satellite apogee and perigee height vs true anomaly at the second cut and 
libration angle at the first cut, for 0 M =59 deg. 
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Apogee Height for Theta Max = 60 deg 



52 




Apogee Height for Theta Max = 61 deg 
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Dynamics and Control of SEDSAT 


Introductory Remarks 

The satellite SEDSAT, built by a team of students of the University of Alabama at 
Huntsville, will be injected into a higher orbit from the Space Shuttle in July 1997 by 
using a 20-km-long librating tether acting as a sling shot. Since the satellite is required to 
have a lifetime of about 3 years, without reboosting, the height and eccentricity of the 

final orbit are a primary concern. The injection velocity (AV) at the satellite release must 
be maximized in order to maximize the apogee height of the final orbit and consequently 
the satellite lifetime. Moreover, the final orbit must be as insensitive as possible to 
variations of the deployer's friction parameters in order to provide an accurate orbital 

transfer of the satellite. Consequently, the magnitude, direction, and timing of the AV at 
satellite release (i.e., when the tether crosses the local vertical and the tether is cut from 
the Shuttle) must be accurate and as insensitive as possible to variations of the model 
parameters. A non-linear control law (feedforward-feedback) for driving the satellite to 
the desired state at the end of deployment was developed during the 1994-1995 research 
activity on this grant [1]. The control law is based on predefined reference profiles which 
are dependent upon the conditions at satellite ejection from the Shuttle. These conditions 
at satellite ejection have changed since the control law that was delivered to NASA in 
December 1994 due to Shuttle safety considerations. Consequently, new reference 
profiles had to be generated and some of the control parameters modified in order to cope 
with the new initial conditions. 

New Reference Profiles and Control Parameters 

Before describing the changes to the SEDSAT deployment control law, we will 
recall the basics of the SEDS deployer tension model and the control law. The reader 
should read Refs. [1-2] for a more detailed treatment of the subject. 

SEDS Deplover's Tension Model 

The tension model for the SEDS deployer is as follows: 


T = 





(92) 


where A re l = 1 - A L/L en d> L enc i = final tether length = 20 km, A = tether annulus solidity 
= 0.9424, E = area exponent = 0.6, B = brake parameter = 2nfn , n = number of brake 
turns,/ = friction coefficient = 0.18 (best estimate), To = static (or minimum) tension, I = 

inertia multiplier = 4.1, P = linear density of tether = 3.3 kg/km. Go = null exit angle, 

and 0 = tether’s exit angle (for a deployer aligned along the Nadir, 0 coincides with the 
in-plane libration angle). In summary, this approximate tension model consists of a static 

tension term To, an hydraulic tension term proportional to Z?, an exponential brake 
multiplier, and an exponential exit-angle multiplier. Parameters of the tension model are 
uncertain as they depend on environmental conditions and many other uncontrollable, or 
poorly controllable, factors. The most influential, from the deployment dynamics 
standpoint, and also less certain parameter of the tension model is To ■ A likely estimate 
of the value of To from the flights of SEDS-I and II is about 20 mN with an expected 
variability of ±10 mN. 
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Basics of Deployment Control Law 


SEDSAT is ejected from the Shuttle with initial conditions (i.e., ejection angle and 
magnitude of ejection velocity) such that the amplitude of the final tether libration in the 
orbital plane is insensitive to variations of the static tension of the deployer's tension 
model. During this time (free-deployment phase), the brake is not engaged. At a tether 
length greater than about 15 km, the brake is engaged to slow down the tether exit 
velocity. The brake utilizes a feedback-feedforward, non-linear control law [2-4] to lead 
the system to the desired state vector at the end of deployment. The control law is based 
on precomputed length and length rate profiles (i.e., the non-linear portion of the control 
law) and on a linear feedback which keeps the system state along the desired state 
trajectory. 

The timing of the crossing of the local vertical is dependent upon the final libration 
amplitude. This amplitude is affected by the errors of the conditions at ejection with 
respect to the desired ejection conditions. A change of the magnitude of the ejection 
velocity require the computation of new reference profiles and also the evaluation of the 
new ejection angle that provides insensitivity to the static tension during the free- 
deployment phase. 

The ejection velocity of SEDSAT has been changed several time from the original 
Vo = 1.47 m/s to the latest Vo — 2.5 m/s as a consequence of Shuttle's safety 
considerations. A number of reference profiles (from SEDSAT_Ref39 to 
SEDSAT_Ref53) have been derived to accommodate these changes. As a result, the 

ejection angle must be changed from 0o = -31° (upward and forward) consistent with Vo 

= 1.47 m/s to Go = - 21° consistent with Vo = 2.5 m/s. 

In this report we present only the results of the latest reference profile and the 
associated control parameters, namely, SEDSAT_Ref53_lSep95 (the numbering is 
sequential from the beginning of the SEDSAT control law analysis). At the time of the 
delivery of the control law to NASA/MSFC, on 1 September 1995, the system and 
tension nominal parameters were as follows: 

Tension model parameters 

To = minimum tension = 20 mN 
I = inertia multiplier = 4.1 
A = annulus solidity = 0.9424 
E = area exponent = 0.6 
f = friction coefficient = 0. 1 8 

System parameters 

Satellite mass = 36.3 kg 
Tether linear density = 0.33 kg/km 
Shuttle orbit = 297 km circular 
Orbital inclination = 57° 

Tether length = 20 km 
Tether diameter = 0.8 mm 
Tether material = Spectra 1000 
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Table 14. Control and system parameters for SEDSAT_Ref39 and SEDSAT_Ref53 


CONTROL PARAMETERS 


No. 

PARAMETER 

TYPE 

Units 

REF39 

REF53 (new) 

1 . 

c 

Filter Coefficient 

(none) 

0.125 

0.125 

2. 

K1 

TumCount Gain 

(1/Tum) 

0.002 

0.002 

3. 

DZTC 

TumCount Deadzone 

(Turn) 

5 

5 

4. 

TCELIM 

Maximum Turn 
Count Error 

(Turn) 

5000 

6500 

5. 

K2 

Turn Count 
Rate Gain 

(s/Tum) 

0.2 

0.2 

6. 

DZTCR 

TumCountRate 

Deadzone 

(Tum/s) 

0.5 

0.5 

7. 

TCRELIM 

Maximum Turn 
CountRate Error 

(Tum/s) 

10 

25 

8. 

WAILP 

Wrap Increment 
Upper Limit 

(none) 

2 

2 

9. 

TBDs 

Time after which 
bias is applied 

(s) 

65535 

65535 

10. 

BIAS 

Brake post Bias 

(Turn) 

0 

0 

11. 

WACLP 

WrapAdjustment 
Upper Limit 

(Turn) 

7 

7 

12.* 

TCBS 

TumCount BrakeStop 

(Turn) 

46000 

46000 

13* 

A1 

Coeff_l in 
Variable Gains 

(none) 

0.6980403 

0.6980403 

14.* 

A2 

Coeff_2 in 
Variable Gain 

(none) 

5.73803 le-6 

5.73803 le-6 

15. 

STOPDEPLOY 

Time from eject, for 
Brake Ramping up 

(s) 

5000 

4600 


— INITIAL CONDITIONS AT EJECTION 

— 



Vo 

Ejection velocity 

(m/s) 

1.47 

2.5 


00 

Ejection angle 
(up and forward) 

(deg) 

-31 

-21 


* These parameters are very likely to change for the SEDSAT flight tether 
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Table 14 shows the control law parameters for SEDSAT_Ref53. The logical names 
of the control variables in the first column of Table 14 are in accordance with the 
document SEDS Deploy er Flight Software Requirement Specification (DM06) [7]. 

Table 14 not only shows the values of the control and system parameters for the latest 
control law SEDAT_Ref53_lSep95 (Vo = 2.5 m/s) but also shows the corresponding 
values of the previously released control law SEDSAT_Ref39_6Dec94 (Vo = 1.47 m/s). 

As indicated in Table 14, some of the control parameters will likely need to be 
changed as the characteristics of the flight tether for SEDSAT will become known after 
the deployment tests are conducted. 

The reference profiles, which are expressed in terms of number of deployed turns 
and turn rate, will also have to be modified as a consequence of a change in the 
relationship between tether length and tether turns. The maximum number of tether turns 
on the spool used to derive SEDSAT_Ref53 is 46200. The reference tables for 
SEDSAT_Ref53 are shown in Appendix I. The reference length, rate and brake profiles 
are shown in Fig. 31. Note that the in-plane libration angle and angular rate, depicted in 
Fig. 31, are not reference variables. They are shown here only for the sake of 
completeness. 

The dynamics of SEDSAT deployment according to the new reference profiles, for 
reference conditions, is shown in Fig. 32. Under reference conditions, the satellite will 
take 4480 s to reach the final tether length of about 20 km and 4780 s to cross the local 
vertical (LV). 
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Fig. 31. Reference profiles for SEDSAT_Ref53. 
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Fig. 32. Deployment dynamics for SEDSAT_Ref53 (reference conditions) 
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Sensitivity of New Deployment Control Law to Errors 

We rewrite here below the reference conditions which were adopted for deriving 
SEDSAT_Ref53: 

Tension model parameters 

To = minimum tension = 20 mN 
I = inertia multiplier = 4. 1 
A = annulus solidity = 0.9424 
E = area exponent = 0.6 
f = friction coefficient = 0.18 

Initial conditions 

Vo = ejection velocity = 2.5 m/s 

00 = ejection angle = -21° (up and forward) 

Table 2 shows the timing errors At at the crossing of LV with respect to the nominal 
value for errors affecting the minimum tension To, the Inertia Multiplier I, the ejection 

angle 6o and the ejection velocity Vo- 

Table 2. At (s) errors at the crossing of LV for SEDSAT_Ref53 


Lead cases 

At error at LV 

At error at LV 

Lag cases 

Ref conditions but 
To = 10 mN (-50%) 

-8 s 

+8 s 

Ref conditions but 
T 0 = 30 mN (+50%) 

Ref conditions but 
1 = 5.1 (+25%) 

-40 s 

+60 s 

Ref conditions but 
1 = 3.1 (-25%) 

Ref conditions but 
0 O = -23° (+10%) 

-54 s 

+60 s 

Ref conditions but 
0 O = -19° (-10%) 

Ref conditions but 
V 0 = 2.63 m/s (+5%) 

-90s 

+94 s 

Ref conditions but 
V 0 = 2.38 m/s (-5%) 


Timing errors due to initial conditions are proportional to initial condition errors as 
follows: 

18 s per 1% of ejection velocity (magnitude) error 

28 s per 1° (or 5%) ejection angle error 

Timing errors at LV crossing (and release) of less than ±100 s have negligible effect 
(< 100 m altitude variations) on the perigee height after release. The effect on the apogee 
altitude is less than ±1% of the apogee height after release. Timing errors of ±100 s 
correspond to maximum tether angle deviations with respect to LV of about ±10°. 

In Figs. 33(a)-33(h) the dynamics of SEDSAT for off-nominal values of the 
parameters above is compared to the nominal deployment dynamics. 
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Fig. 33(a). Sensitivity of deployment to Minimum Tension 
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Fig. 33(b). Sensitivity of deployment to Minimum Tension (total AV) 
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Fig. 33(c). Sensitivity of deployment to Inertia Multiplier 
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Fig. 33(e). Sensitivity of deployment to Ejection Angle 
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Fig. 33(f). Sensitivity of deployment to Ejection Angie (total AV) 
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Fig. 33(g). Sensitivity of deployment to Ejection Velocity 
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Fig. 34. SEDSAT deployment for Vo = 2 m/s (ejection mechanism ma lfunction). 
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Fig. 34(b). SEDSAT deployment for Vo = 2 m/s (malfunction; total AV). 
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Finally, a malfunction of the ejection mechanism was analyzed. In this scenario, the 
ejection velocity is reduced from Vo = 2.5 m/s to Vo = 2 m/s due to an assumed breakage 
of one ejection spring. The dynamics of SEDSAT for this malfunction is shown in Figs. 
34(a)-34(b). The satellite deployment still reaches a length of about 20 km (19.8 km to 
be specific) in 4488 s (only 8 s longer than the nominal case) but the tether will cross LV 
at t = 5182 s (i.e., 402 s longer than the nominal case). The orbital injection of the 
satellite would still be acceptable if the satellite is released at the time the tether actually 
crosses LV. 

The AV imparted to the center of mass of the tether-satellite system consists of a 
hanging component (that depends on tether length) and a swinging component that is due 
to the motion of the tether-satellite system with respect to the LH-LV reference frame. 
The hanging component is less affected by the timing errors than the swinging 
component. The (nominal) hanging component AVh is equal to 21 m/s at release while 
the (nominal) swinging component AV$ is equal to 32 m/s. The nominal (i.e., reference 

conditions) post-release orbit of the satellite-tether system is about 317 km X 560 km, 
with respect to a spherical Earth of equatorial radius, after the first tether cut at the 
Shuttle end. 

The absolute time of the first cut (or equivalently the orbital anomaly) is not critical 
for SEDSAT because there is no specific requirement on the phase of the final orbit for 
this mission. However, the sensitivity of the release time to variations of the deployer's 
model parameters can be more or less important depending on the level of crew and on- 
board equipment involvement in estimating the tether LV crossing. In the simplest 
scenario, the LV crossing could be based solely on the deployment start time and, 
consequently, the insensitivity of the maneuver duration to different deployer's 

parameters is important for obtaining an accurate AV at release. 

For example the timing errors associated with the expected variations of the static 

tension is only ±10 s. The hydraulic component (proportional to Z?) in the deployer's 
tension model also has an effect on the deployment dynamics and, hence, on the timing 

and magnitude of AV at the crossing of LV. Unfortunately, the ejection angle that 
cancels the effect of the hydraulic tension term has the opposite sign of the ejection angle 
that cancels the effect of the minimum tension. It can be inferred from Fig. 33(c), which 
show the libration angle for different values of the hydraulic tension term, that a variation 
of the inertia multiplier I of ±25%, with respect to the reference value of 4.1, implies a 
variation of the timing to cross LV of about ±60 s. 

The variation of the horizontal component of the total AV at release due to timing 
errors of ±100 s is less than ±1 m/s at the nominal time of LV crossing . The variation of 

the vertical component of the total AV is less than ±7 m/s at the nominal time of LV 
crossing . However, a vertical AV produces an increase (or decrease) of the orbital height 
which is about 1/4 of the effect produced by an horizontal AV. 

In summary, if the tether cut is done at the nominal time of LV crossing from the 
start of deployment, the dispersion of the horizontal AV component is minimized If the 

same tether cut is done at the actual crossing of LV, the dispersion of the vertical AV 
component is minimized . The latter case requires more crew and on-board equipment 
involvement. The two cases produce comparable dispersion in terms of the height of the 
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final orbit. The former case can produce a slightly different phasing of the initial and 
final orbit which, however, is not an issue for SEDSAT at present. 

The time accuracy at release is definitely important for future missions which are 
expected to use a SEDS-type deployer for the atmospheric reentry of capsules. In this 
case the dispersion of the absolute time at release determines to a large extent the 

dispersion of the reentry capsule footprint. For instance, a timing error of the reentry AV 
in LEO of ±20 s implies, approximately, a ±148 km (±80 nmi) error in the length of the 
footprint. It is worth remembering that the length of the footprint of the Gemini capsule 
was typically ±185 km (±100 nmi) [7]. 

Concluding Remarks 

In 1997, SEDSAT will the deployed on a 20-km-long tether from the Shuttle orbiting 
on a 297 -km circular orbit and released into a higher orbit. The final tether length of 20 
km is reached in 1 hr 14 min 40 s according to the new control law SEDSAT_Ref53 
described in this report. Upon crossing the local vertical at t = 1 hr 19 min 40 s (i.e., 
5 min after reaching the final deployed length), the tether is cut at the Shuttle and, 
consequently, the satellite with the attached tether is injected into an higher orbit thanks 
to the librating tether acting as a sling shot. After the first tether cut (at the Shuttle end) 

the orbit of the satellite-tether system will be about 317 km X 560 km, with respect to a 
spherical Earth. 

The expected timing errors at LV crossing (and release) of less than ±100 s have 
negligible effect (< 100 m altitude variations) on the perigee height after release. The 
effect on the apogee altitude is approximately ±1% of the apogee height after release (i.e., 
about 6 km). Timing errors of ±100 s correspond to maximum tether angle deviations 
with respect to LV of ±1 1°. 
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Appendix A: RITZTRIG.MS 


RITZTRIG.MS 

RITZ METHOD WITH TRIGONOMETRIC FUNCTIONS 

> b:=1/(1 +32/(3.3E-4*20E3))/2; 

> with(linalg): 

b := .08549222800 

Warning: new definition for norm 

Warning: new definition for trace 


> F:=v->K*cos(2*n*Pi*v)+sin(2*n*Pi*v); 

F :=v ->K cos( 2 n n v ) + sin( 2 n n v ) 

> R2:=int(F(v),v=b-1..b)+(1/2/b-1 )*F(b): 

> Kn2:=solve(R2=0,K) $ n=1..4; 

Kn2 := -.5955804727, -1.845948652, 24.55995647, 1.533481529 

> DF2:=diff(diff(F(v),v),v): 

> R1:=subs(v=b,DF2); 

> Kn1:=solve(R1=0,K) $ n=1..4; 

R1 := -4 K cos(. 1709844560 n n) n 2 n 2 -4 sin(. 1709844560 nn) n 2 it 2 

Knl := -.5955804728, -1.845948652, 24.55995647, 1.533481529 

> evalf(subs({K=Kn1[1],n=1},R1)); 

> evalf(subs({K=Kn1 [1 ],n=1 },R2)); 

> evalf(subs({K=Kn2[1],n=1},R1)); 

> evalf(subs({K=Kn2[1 ],n=1 },R2)); 

0 

0 

0 

0 

> evalf(subs({K=Kn1[2],n=2},R1)); 

> evalf(subs({K=Kn1[2],n=2},R2)); 

> evalf(subs({K=Kn2[2],n=2},R1)); 

> evalf(subs({K=Kn2[2],n=2},R2)); 

0 

-.1 10‘ 8 

0 

-.1 10" 8 

> evalf(subs({K=Kn1[3],n=3},R1)); 

> evalf(subs({K=Kn1[3],n=3},R2)); 

> evalf(subs({K=Kn2[3],n=3},R1)); 

> evalf(subs({K=Kn2[3],n=3},R2)); 

0 

.1 10' 8 

0 

.1 10' 8 
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> evalf(subs({K=Kn1[4],n=4},R1)); 

> evalf(subs({K=Kn1 [4],n=4},R2)); 

> evalf(subs({K=Kn2[4],n=4},R1)); 

> evalf(subs({K=Kn2[4],n=4},R2)); 

.1 10' 6 
-.1 10' 8 
.1 10‘ 6 
-.1 10' 8 

> FI :=subs({K=Kn1 [1],n=1 },F(v)); ' 

FI := -.5955804728 cos( 2 it v ) + sin( 2 it v ) 

> F2:=subs({K=Kn1[2],n=2},F(v)); 

F2 := -1.845948652 cos(4 jt v ) + sin(4 it v ) 

> F3:=subs({K=Kn1[3],n=3},F(v)); ~ 

F3 := 24.55995647 cos(6 it v) + sin(6 it v) 

> F4:=subs({K=Kn1[4],n=4},F(v)); “ 

F4 := 1.533481529 cos(8 it v) + sin(8 it v) 

> plot({F1 ,F2,F3,F4},v=b-1 ..b); 



>c:=array(1..4,1..4); 

c :=array(l .. 4, 1 ..4, [ ]) 

> c[1,1]:=int(F1*F1,v=b-1..b); " 

> c[1 ,2]:=int(F1 *F2,v=b-1 ..b);c[1 ,2]:=0; 

> c[1 ,3]:=int(F1 *F3,v=b-1 ,.b);c[1 ,3]:=0; 

> c[1 ,4]:=int(F1 *F4,v=b-1 ..b);c[1 ,4]:=0; 

c -=.6773580491 

A 


c := .3395305451 10' 9 

1,2 
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c := -.5570423009 10' 8 

1,3 



c := .2232413334 10‘ 9 

1,4 


> c[2,1]:=int(F2*F1,v=b-1..b);c[2,1]:=0; 

> c[2,2]:=int(F2*F2,v=b-1 ..b); 

> c[2,3]:=int(F2*F3,v=b-1..b);c[2,3]:=0; 

> c[2,4]:=int( F2*F4 , v=b- 1 ..b);c[2,4]:=0; 




.3395305451 10' 9 



:= 0 


c := 2.203763211 

2,2 


2, 3 " 
2, 3 :i 


c -=-.5305164769 10‘ 9 

2,4 


> c[3,1]:=int(F3*F1 ,v=b-1 ..b);c[3,1]:=0; 

> c[3,2]:=int(F3*F2,v=b-1 ..b);c[3,2]:=0; 

> c[3,3]:=int(F3*F3,v=b-1 ..b); 

> c[3,4]:=int(F3*F4,v=b-1 ..b);c[3,4]:=0; 




:= 0 


-.5570423009 10‘ 8 




c := 302.0957310 
3,3 


c := -.1 136821022 10 -8 

3,4 


> c[4,1]:=int(F4*F1,v=b-1..b);c[4,1]:=0; 

> cf4,21:=int(F4 > F2,v=b-1 ..b);cr4,21:=0; 
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> c[4,3]:=int(F4*F3,v=b-1..b);c[4,3]:=0; 

> c[4,4]:=int(F4*F4,v=b-1..b); 

c : = .2232413334 10' 9 
4, 1 


c „ .--.5305164769 10" 9 
4,2 


c „ :=-. 1136821022 10“ 8 
4, 3 


c r= 1.675782800 
4,4 


> cp:=array(1 ..4,1 ..4); 

cp := array( 1 „ 4, 1 .. 4, [ ]) 

> cp[1,1J:=int(Frdin(Fi,v)’v,v=b-l..b); 

> cpfl ,2]:=int(Frdiff(F2,v)*v t v=b-1 ..b); 

> cp[1 ,3]:=int(F1*diff(F3,v)*v,v=b-1 ..b); 

> cp[1 ,4]:=int(Frdiff(F4,v)*v,v=b-1 ..b); 


cp := -.3386790252 

I , I 

cp, ^ := 1.629034753 
1,2 

cp -=-10.72859139 

I, J 

cp, .= -.5682202935 
1,4 

> cp[2,1]:=int(F2*diff(F1 ,v)*v,v=b-1 ..b) 

> cp[2,2]:=int(F2*diff(F2,v)'v,v=b-1 ..b) 

> cp[2,3]:=int(F2*diff(F3,v)*v,v=b-1 ..b) 

> cp[2,4]:=int(F2*diff(F4,v)*v,v=b-1 ..b) 

° P 2 1 := -1.629034753 

cp := -1.101881606 

2, 2 

cp „ := -61.92499786 
2, 3 

cp 2 4 := -2.562300523 

> cp[3,1]:=int(F3*diff(F1 ,v)*v,v=b-1 ..b) 

> cp[3,2]:=int(F3*diff(F2,v)*v,v=b-1 ..b) 

> cp[3,3]:=int(F3*diff(F3,v)*v,v=b-1 ..b) 

> cp[3,4]:=int(F3*diff(F4,v)*v,v=b-1 ..b) 
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cd := 10.72859140 

3, 1 

CD „ := 61.92499789 
3, 2 

cp„ := -151.0478655 

3, 3 

CD := 77.14261558 

> cp[4, 1 ]:=int(F4*diff(F1 ,v)*v,v=b-1 ..b); 

> cp[4,2]:=int(F4*diff(F2,v)*v,v=b-1 ..b); 

> cp[4,3]:=int(F4*diff(F3,v)*v,v=b-1 ..b); 

> cp[4,4]:=int(F4*diff(F4,v)*v,v=b-1 ..b); 

CD , := .5682202936 

4, 1 

CD := 2.562300523 

*5 ^ 

cp „ := -77.14261558 
4,3 

cp -=-.8378914009 

4>4 

> cs:=array(1 ..4, 1 ..4): ~ 


> cs[1 ,1 ]:=int(F1 *diff(diff(F1 ,v),v)*((b-1 )**2-v**2),v=b-1 ..b); 

> cs[1,2]:=int(Frdiff(diff(F2,v),v)*((b-1)**2-v**2),v=b-1..b); 

> cs[1 ,3]:=int(F1 *diff(diff(F3,v),v)*((b-1)**2-v**2),v=b-1 ..b); 

> cs[1,4]:=int(Frdiff(diff(F4,v),v)*((b-1)**2-v**2),v=b-1..b); 

cs := -15.87987864 

1 , X 

cs := 8.688185347 

1,2 

cs, „ := -48.27866129 

1.3 

cs, := -2.424406593 

1.4 

> cs[2,1]:=int(F2*diff(diff(F1,v),v)*((b-1 )"2-v**2),v=b-1 ..b); 

> cs[2,2]:=int(F2*diff(diff(F2,v),v)*((b-1)**2-v**2),v=b-1 „b); 

> cs[2,3]:=int(F2*diff(diff(F3,v),v)*((b-1)**2-v**2),v=b-1 „b); 

> cs[2,4]:=int(F2*diff(diff(F4,v),v)*((b-1)**2-v*’2),v=b-1..b); 

cs := 2.172046338 

2 , 1 

cs := -203.3531074 

2,2 

cs „ := -445.8599847 
2,3 

cs _ := -13.66560282 
2, 4 

> cs[3,1]:=int(F3*diff(diff(F1,v),v)*((b-1)**2-v**2),v=b-1..b); 

> cs[3,2]:=int(F3*diff(diff(F2,v),v)*((b-1 )**2-v**2),v=b-1 ..b); 

> cs[3,3]:=int(F3*diff(diff(F3,v),v)*((b-1 )**2-v**2),v=b-1 ,.b); 


A 5 


Appendix A: RITZTRIG.MS 


> cs[3,4]:=int(F3*diff(diff(F4,v),v)*((b-1)**2-v**2),v=b-1..b); 

cs„ , := -5.364295710 
3, 1 

cs := -198.1599933 

L 

cs -=-62532.19241 

3,3 

cs„ , := 705.3039141 

3, 4 

> cs[4,1]:=int(F4*diff(diff(F1,v),v)*((b-1)**2-v**2),v=b-1..b); 

> cs[4,2]:=int(F4*diff(diff(F2,v),v)*((b-1)**2-v**2),v=b-1 ..b); 

> cs[4,3]:=int(F4*diff(diff(F3,v),v)*((b-1)**2-v**2),v=b-1 ..b); 

> cs[4,4]:=int(F4*diff(diff(F4,v),v)*((b-1 )‘*2-v**2),v=b-1 ..b); 

cs, := -.1515254114 

4, 1 

cs, := -3.416400702 

4.2 

cs, „ := 396.7334516 

4.3 

cs, ,:= -616.0203632 

4.4 

> ai:=array([a1,a2,a3,a4]); ~ ‘ 

ai := [al a2 a3 a4] 

> ais:=array([a1s,a2s,a3s,a4s]); ~ ~ ~~ 

ais := [als a2s a3s a4s ] 

> A:=evalm(inverse(c)&*(cs/2-cp)); 

" -11.22192363 4.008305393 -19.79859733 -.9507276165* 

A = 1.232009822 -45.63769444 -73.05911707 -1.937822023 

-.04439234943 -.5329601775 -102.9973122 .9119934944 

.-.3842878678 -2.548361800 164.4063546 -183.3007775 
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RITZLEGO.MS 

n omin ation of the Constants Ki in trial solution not have Legendre 

> b:=1/(1+32/(3.3E-4*20E3))/2; 

> epsilon:=b/(1-b); 

b := .08549222800 
e := .09348441929 


> with(orthopoly): 

Maple does not have Legendre functions of the second kind, but they can be given in 
terms of Pn(x) through the formula (cf. Abramowitz, 8.6.19): 


> WW:=(n,x)->sum(P(i-1,x)*P(n+1-i,x)/i,i=1..n+1); 

71+1 


ww 


XT' P(i — 1,*) P(tz+ 1 -i,x) 
: =( n '*)^ 2, i 


i = 1 


> Q:=(n,x)->P(n,x)*ln((x+1 )/(1-x))/2-WW(n-1 ,x); 

1 ( x + O 
Q:=(n,x)->-P(7t,x)ln - — - 

2 1 — X 


-WW(n-U) 


lambda0:=n->n*(n+1)/2; 


A.0 :=ti ->-ti (ti + 1) 


> lambdal :=n->4*(1 -lambdaO(n))’(n+1/2)*GAMMA(1+n/2)**2/GAMMA((1+n)/2)**2/Pi/lambdaO(n); 

l\ ' ' 2 

(1 - X0(n)) I 7t +~ 


X»1 ti — > 4 ■ 


-y 


r ( i+ i n 


i 


\ »+^J JCAjO(7 l) 


The solution of the DifLEq. in space is: 

> F0:=(n,x)->P(n,x); 

> F1:=(n,x)->P(n,x)*(K+2*lambda1(n)*int(P(n,u)*Q(n,u),u=0..x))+ 

> Q(n,x)*((lambdaO(n)-1)*4*GAMMA(1+n/2)**2/lambdaO(n)/Pi/GAMMA((1+n)/2)**2 

> -2*lambda1 (n)*int(P(n,u)**2,u=0..x)); 

> F:=F0+epsilon*F1; 

F0 := P 


FI := (n,x) -+P(ti,x) 


K+2\\(n) ^ P(ti, u) Q(n,u) du 


V 


+ Q(ti, x) 


(X0(7i)-i)rfi+^7tj 


J 1 1> 

X0(n) n I"l 2 71 + 2 


0 

2 


— 2 A,l( 


n) P (ti, u 

0 


) 2 du 


F:=P+ .09348441929 FI 
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Where the constant K Is not determinedby the B.C.; probably we would have to go to the 
e*e approximation to have it determined. We will try later to choose a value that minimizes 
deviation from the approximate B.C. 

> Iambda:=n->lambda0(n)+epsilon*lambda1 (n); 

X :=n ->X0(n) + eXl(n) 

The Diff.Equation is: ~ ~ 

> ED:=(n,x)->(1-x**2)*diff(diff(F(n,x),x),x)-2*x*diff(F(n,x),x)+2*lambda(n)*F(n,x); 

ED := 

(n,x) — >(l -3c^)diff(diff(F(n > j:) > 3i:) > 3c)-2xcliflF(F( n, x), x) + 2 X(n) F (n,x) 

We will see how close the asymptotic development solution meets the equation; since K 
is not determined, we will give several values: 

> plot({subs(K=0, ED(5,x)),subs(K=0.5,ED(5,x)),subs(K=1,ED(5,x))},x=-1.. epsilon, title=' Residue for n=5') 
> ; 


Residue for n=5 



> plot({subs(K=0,ED(3,x)),subs(K=0.5,ED(3,x)),subs(K=1 ,ED(3,x))},x=-1 ..epsilon, title=" Residue for n=3' ) 



> plot({subs(K=0,ED(1 ,x)),subs(K=0.5,ED(1 ,x)),subs(K=1 ,ED(1 ,x))},x=-l ..epsilon,title=' Residue for n=1') 

> ; 
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We will check that the solutions FO and FI are of the same order: 

> plot({F0(5,x),subs(K=0, FI (5, x)),subs(K=0.5, FI (5, x)),subs(K=1,F1 (5, x))},x=-1.. epsilon, title=' Functions f 

> or n=5 and different K'); 

Functions for n=5 and different K 



>plot({F0(3,x),subs(K=0,F1(3,x)),subs(K=0.5,F1(3,x)),subs(K=1,F1(3,x))},x=-1.. epsilon, title=' Functions f 
> or n=3 and different K'); 
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Functions for n=3 and different K 



> plot({FO(1 ,x),subs(K=0,F1 (1 ,x)),subs(K=0.5,F1 (1 ,x)),subs(K=1 ,F1(1 ,x))},x=-1 ..epsilon, title='Functions for n= 

> 1 and different 1C); 

Functions for n=l and different K 



The three B.C. (second derivative zero at b/(l-b), first derivative=lambda*function/b, 
and Center of Mass integral) are the same when the Diff.Eq. is fulfilled; nnw this 
solution is approximate, the three B.C. will also be approximately met. We will try to 
choose the value of K so as to bring them closer. 

> phi1:=eval(F(1,x)); 

<|>1 := x + .09348441929 x K 

> phi3:=eval(F(3,x)): 

> phi5:=eval(F(5,x)): ~~ “ — ■■ 


> phi1s:=diff(diff(phi1,x),x); 

phils := 0 

> phi3s:=diff(diff(phi3,x),x): 

> phi5s:=diff(diff(phi5,x),x): 

> phil p:=drff(phi 1 ,x); 
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philp := 1 + .09348441929 K 

> ph i3p: =d itt( piTrJTxJ; 

> phi5p:=ditt(phib,x): 

Residuals for tbe different B.C.: 

> R1s:=subs(x=epsilon,phi1s); 

> R3s:=evalf(subs(x=epsilon,phi3s)); 

> R5s:=evalf(subs(x=epsilon,phi5s)); 

Rls ~ 0 


R3s := -.1325489892 + .1310900497 K 
R5s := .1697377198 - .4467859534 £ 

> R1p:=evalf(subs(x=epsilon, phi1p’epsilon-(lambda(1))’phi1)); 

> R3p:=evalf(subs(x=epsilon, phi3p*epsilon-(lambda(3))*phi3)); 

> R5p:=evalf(subs(x=epsilon, phi5p*epsilon-(lambda(5))*phi5)); 

Rip := 0 


R3p := .00209883271 + .05704708175 K 
R5p := -.1549035234 - .1949139959 

> R1i:=evalf((1-b)~int(F(1,x),x=-1..epsilon)+(1/2/b-1)’F(1, epsilon)); 

> R3i:=evalf((1-b)*int(F(3,x),x=-1..epsi!on)+(1/2/b-1)*F(3, epsilon)); 

> R5i:=evalf((1-b)*int(F(5,x),x=-1..epsilon)+(1/2/b-1)*F(5, epsilon)); 

Rli := -.2 10' 9 - .2 10‘ 10 AT 

R3i := -.00433781797 - .1097409326 10‘ 10 /- .05250279174 K 
R5i := .0589231983 + .2194818653 10‘ 10 /+ .07157678227 K 

The first eigenfunction is linear; the value of K does not affect the B.C. at all: make it zero. 

> Kp[1]:=solve(R1p,K); 

> Ks[1]:=solve(R1s,K); 

> Ki[1]:=fsolve(Re(R1i),K); 

Kp x :=K 
Ks ^ :=K 

Ki { := -10.00000000 

> plot(subs(K=o,phii),x=-l.. epsilon, titie= First Eigenfunction tor K=0 ); 
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First Eigenfunction for K=0 



> Kp[3]:=solve(R3p,K); 

> Ks[3]:=solve(R3s,K); 

> Ki[3]:=fsolve(Re(R3i),K); 

> plot({R3p,R3s,Re(R3i)},K=-10..10,title='B.C. Residuals for different 1C); 

B.C. Residuals for different K 



Kp := -.03679123709 
Ks := 1.011129292 
Ki := -.08262071075 

We will choose the one minimizing the residual of the two B.C. which apply to the general 
problem, the integral and second derivative ones: 

> Kav[3]:=fsolve(R3s+Re(R3i),K); 

Kav^ := 1.741844807 

> plot({subs(K=0 1 phi3),subs(K=Kav[3],phi3),subs(K=Kp[3],phi3),subs(K=Ks[3],phi3),subs(K=Ki[3] I phi3)}, 

> x=-1.. epsilon, title='Second eigenfunction for different 1C); 
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> Kp[5]:=solve(R5p,K); 

> Ks[5]:=solve(R5s,K); 

> Ki[5]:=fsolve(Re(R5i),K); 

> plot({R5p,R5s,Re(R5i)} I K=-10..10,title='B.C. Residuals for different FC); 

B.C. Residuals for different K 



Ki 5 := -.8232166414 

> Kav[5]:=fsolve(R5s+Re(R5i),K); 

Kav := .6094225187 
5 

> plot({subs(K=0,phi5),subs(K=Kav[5],phi5),subs(K=Kp[5],phi5),subs(K=Ks[5],phi5),subs(K=Ki[5],phi5)}, 

> x=-1.. epsilon, title='Third Eigenfunction for different K'); 
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Third Eigenfunction for different K 



The value of K does not affect much the shape of the 3rd and 5th Eigenfunctions; the 
value of the B.C. residuals is small, except for the 5th function. So we will choose K=0 
for the 1st Eigenfunction, and K=Kav for the 3rd and 5th, which minimizes the sum of 
the two ideally identical B.C. If we try to minimize the first derivative B.C., the function 
will be closer to the eigenfunctions, but the B.C. Ritz trial functions must fulfill would be 
less accurately met; so we choose K to minimise only Ri and Rs. 

> plot({subs(K=0, phi 1),subs(K=Kav[3],phi3),subs(K=Kav[5],phi5)},x=-1..epsilon,title=' Functions for K min 

> imizing B.C. residuals'); 

Functions for K minimi zing B.C. residuals 



> phi1:=subs(K=0,phi1); 


<J>1 := X 


> phi3:=subs(K=Kav[3], phi3); 


$3 := | x 3 - 1 x + .09348441929 ^|x 3 xj ^ 1.741844807 ln^- 


505 4 
— -x 
32 


64 

32 \ -1 + x 


5 45 2 15 1 

x J + — x*——]n(2) +.09348441929 

O 8 J 
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<(>3 := ^ x 3 — ^ x + .09348441929 f|x 3 -|x 
Y 2 2 1 2 2 , 


1.741844807 


315 , ( 1 + x 

64 X -1 + x ) 


16 ^ -1+x 


505 4 
32 * 


3 375 6 375 

+ 32 * 64 



)* 7 +^lr 

f 2 1 ] 

1 -1+x, 

r + 8 k 

1 -1+x) 


315 f 1 +x A 
+ -rr- ln| - 


32 


-1 + x 


5 45 2 15. ... 

x +yx -yln(2) 


+ .0934844 1929 f^f^x 3 -^x 
2 V 2 2 . 


1+x) 5 2 2 N 

1 — x J 2 X + 3 . 


15 375 7 315 3 315 5 

8 + 32 * + 32 X 16 X 


\ 

> phi5:=subs(K=Kav[5],phi5); 

63 5 35 3 15 f 63 5 35 3 15 V 

:= — x - — x " 3 + — x + .09348441929 — x - — x J + — x .6094225187 

8 4 8 l 8 4 8 A 


♦ 5: =8 
105 
32 


{ 2 1 ] 

992985 8 

86625 

{ 1+ *1 

1 -1+*J 

1 2048 X 

“ 4096 k 

l -1+x) 


l+i) 3 278495 4 

I x x 


2048 


410/4D 1+X 

4096 \ -1+x 

121275 , ( 1+x 


1024 

105 
32 


■ In 


-1 +x ) 


105 , ( 1+x 

——In 

64 -1+x 

560175 , ( 1+x 


x 1 1 - Arr ln| 

J 

' 5 

* 2048 


1024 X - 


1 + x 


Id - 


f 


ln(2) +.09348441929 


-1 + x 

35 
4 


^7 525 2 416745 10 826399 6 
x + — x^+ x iV3 + x° 


32 


2048 


2048 


1 f63 5 35 3 15 'l ( 1+x) 21 4 9 2 9 

l{T X ~T X + T x r[jxj-T x + 2 X -20 


3 f 

r S 3 

3 ) 

1 | 

r 3 2 

n 

2 

— x 




-x - 

— 


4 1 

,2 

2 J 

3i 

,2 


/ 


f 105 121275 5 

t 32 512 X + 2048 * 


-- -x 2 -- 

3 v2 2 )) 

86625 3 560175 7 282975 9 416745 n') 

2048 * + 1024 X 512 X + 2048 * J 
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RITZLEG1.MS: Ritz matrix coefficients 

> assume(x.real) 

> assume(x,real); 

> b:=1/(1+32/(3.3E-4*20E3))/2: 

> epsilon:=b/(1-b): 

> with(orthopoly): 

> WW:=(n,x)->sum(P(i-1 ,x)*P(n+1 -i,x)/i,i=1 ..n+1 ): 

> Q:=(n,x)->P(n,x)*ln((x+1 )/(1 -x))/2-WW(n-1 ,x): 

> lambda0:=n->n*(n+1)/2: 

> lambdal :=n->4*(1 -lambdaO(n))*(n+1/2)*GAMMA(1+n/2)**2/GAMMA((1+n)/2)**2/Pi/lambdaO(n): 

> FO:=(n,x)->P(n,x): 

> FI :=(n,x)->P(n,x)*(K+2*lambda1 (n)*int(P(n,u)*Q(n,u),u=0..x))+ 

> Q(n,x)*((lambda0(n)-1 )*4*GAMMA(1+n/2)**2/lambdaO(n)/Pi/GAMMA((1 +n)/2)**2 

> -2*lambda1(n)*int(P(n,u)**2,u=0..x)): 

> F:=F0+epsilon*F1 : 

> Iambda:=n->lambda0(n)+epsilon*lambda1 (n): 

> Phi1:=F(1,x): 

> Phi3:=F(3,x): 

> Phi5:=F(5,x): 

> phil :=subs(K=0,Phi1 ): 

> phi3:=simplify((subs(K=1 .741 844807, Phi3))): 

> phi5:=simplify(eval((subs(K=0.60942251 87,F(5,x))))): 

> phi1p:=simplify(diff(phi1,x)): 

> phi3p:=simplify(diff(phi3,x)): 

> phi5p:=simplify(diff(phi5,x)): 

> phi1s:=simplify(diff(phi1p,x)); 

> phi3s:=simplify(diff(phi3p I x)): 

> phi5s:=simpirfy(drff(phi5p,x)): 

> with(linalg): 

> A:=matrix(3,3); 

> A[1 ,1]:=Re(int(phi1 *phi1 ,x=-1 ..epsilon)); 

> A[1 ,2]:=Re(int(phi1*phi3,x=-1 ..epsilon)); 

> A[1,3]:=Re(int(phi1*phi5,x=-1.. epsilon)); 

> A[2,1]:=Re(int(phi3*phi1 ,x=-1 ..epsilon)); 

> A[2,2]:=Re(int(phi3*phi3,x=-1 ..epsilon)); 

> A[2,3]:=Re(int(phi3*phi5,x=-1 ..epsilon)); 

> A[3,1]:=Re(int(phi5*phi1 ,x=-1 ..epsilon)); 

> A[3,2]:=Re(int(phi5*phi3,x=-1 ..epsilon)); 

> A[3,3]:=Re(int(phi5*phi5,x=-1 ..epsilon)); 

A := array( 1 .. 3, 1 .. 3, [ ]) 

A := .3336056639 

*■ 9 * 

A :=. 01 171 152560 

1.2 

A := -.005983207937 
1.3 

A , := .01 17 1152560 
2, 1 

A := .2226242651 

2,2 
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A „ „ := -.0006923039365 
2,3 

A „ --.005983207937 

3, 1 

„ := -.0006923039365 
3,2 

,4 -.1260037526 

3.3 

> Ap:=matrix(3,3); 

> Ap[1,1]:=evalf(int(Re(x*phi1*phi1p),x=-1.. epsilon)); 

> Ap[1,2]:=evalf(int(Re(x*phi1*phi3p),x=-1.. epsilon)); 

> Ap[1 ,3]:=evalf(int(Re(x*phi1 *phi5p),x=-1 ..epsilon)); 

> Ap[2,lj:=evalf(int(Re(x*phi3*phi1p),x=-1.. epsilon)); 

> Ap[2,2]:=evalf(int(Re(x*phi3*phi3p),x=-1.. epsilon)); 

> Ap[2,3]:=evalf(int(Re(x*phi3*phi5p),x=-1 ..epsilon)); 

> Ap[3,lj:=evalf(int(Re(x*phi5*phi1p),x=-1.. epsilon)); 

> Ap[3,2]:=evalf(int(Re(x*phi5*phi3p),x=-1 ..epsilon)); 

> Ap[3,3]:=evalf(int(Re(x*phi5*phi5p),x=-1 ..epsilon)); 

Ap := array( 1 .. 3, 1 .. 3, [ ]) 

Ap := .3336056639 

1 , 1 

Ap, := 1.163546054 

1,2 

Ap „ := 1.096924487 

1.3 

Ap := .01 171 152560 
2, 1 

Ap „ := .5937740239 

2,2 

Ap „ := 1.280153230 

2.3 

:= -.005983207929 

3.1 

:= .008282910645 

3.2 

„ := .5252840045 

3.3 

> As:=matnx(3,d); 

> As[1,1]:=evalf(int(Re((1-x*x)*phi1*phi1s),x=-1.. epsilon)); 

> As[1,2]:=evalf(int(Re((1-x*x)*phi1*phi3s),x=-1.. epsilon)); 

> As[1,3]:=evalf(int(Re((1-x*x)*phi1*phi5s),x=-1. epsilon)); 

> As[2,1]:=evalf(int(Re((1-x*x)*phi3*phi1s),x=-1.. epsilon)); 

> As[2,2]:=evalf(int(Re((1 -x*x)*phi3*phi3s),x=-1 ..epsilon)); 

> As[2,3j:=evalf(int(Re((1-x*x)*phi3*phi5s),x=-1.. epsilon)); 

> As[3,1]:=evalf(int(Re((1-x*x)*phi5*phi1s),x=-1.. epsilon)); 

> As[3,2]:=evalf(int(Re((1-x*x)*phi5*phi3s),x=-1.. epsilon)); 

> As[3,3]:=evalf(int(Re((1-x*x)*phi5*phi5s),x=-1 ..epsilon)); 
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As := array( 1 .. 3, 1 .. 3, [ ]) 

-\i :=0 

As, := 2.186056968 

1,2 

As := 2.375025574 

As 2 , 1 ' = ° 

As := -1.267162894 
z 

, := 2.540174068 
2,3 

^ 3 , 1 : =0 

As -.04185944138 

J, Z 

As := -2.380685865 

J) 


J:=evalm(inverse(A)& # (As/z-Ap)); 



- 1.000000000 

-.01663231022 

.0304933737 


-.81 10' 11 

-5.512034854 

-.08915850384 


-.6 10' 10 

.06929390492 

-13.61472320 
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1STCUT.MS 

SEMIAXIS AND EXCENTRICITY PERTURBATIONS: DEPENDENCE ON THETA AND DELTA 
> with(plots): 


DELI A CANNOT BE GREATER THAN THETA 

> delta:=alpha*theta; 

8 := a 0 

> u:=sqrt(3*(sin(theta) A 2-sin(delta) A 2)); 

u := J 3 sin(0)^ - 3 sin(a 0)^ 
PERTURBATTONXTFSEM1AXIS ~ 

> da:=2‘cos(delta)*(u+2); 


da := 2 cos(a 0) / 3 sin(0)^ - 3 sin(a 0)^ + 2 


> plot3d(da,alpha=-1 ..1 ,theta=-Pi/3.“Pi/3); 



EXCENTRICITY: 

> e:=sqrt(3*cos(delta) A 2*(u+1 )*(u+3)+u*u): 

> piot3d(e,aipha=-l..l,theta=-Pi/3..Pi/3); 







MAXIMIZES EXCENTRICITY. LET US SEE THE EFFECT ON PERIGEE HEIGHT: 


> ph:=da-e: 
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AS EXPECTED, THE MAXIMUM PERIGEE DEPENDS ON VERTICAL DISTANCE, AND IS NOT 
AFFECTED BY SPEED. IT IS HIGH AS LONG AS DELTA IS CLOSE TO 
ZERO. APOGEE HEIGHT, HOWEVER, DOES DEPEND ON SPEED: 


> ah=da+e: 
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DECAY1.MS 
DECAY MINIMIZATION 
SHUTTLE IN CIRCULAR ORBIT 
ELLIPTIC ORBIT OF THE CENTER OF MASS 
> with(linalg): 



(epsilon: C.O.M length/circular orbit radius) 

> r:=array([taylor( r0*(1+epsilon*cos(delta)), epsilon), 0,taylor(rO*epsilon*sin(delta), epsilon)]); 

r := [r0~ + r0~ cos(5~) e 0 r0~ sin( 8~ ) s ] 


I ICI 9 1 M IIMMO PiCi 'J M M ljf.1l (Mil IM H Itf 


> v:=array([taylor(-omega0*r0*epsilon*(1+u)*sin(delta),epsilon),0, 

> tay Ior(omega0*r0*( 1 +epsilon*(1 +u)*cos(delta)), epsilon)]); 

v := [ - ©0 r0~ ( 1 + «) sin(6~) 8 0 ©0 + ©0 r0~ ( 1 + w)cos(S~)s] 


i mcv uiti I'Jdi lici h me, rm 


> h:=crossprod(r,v); 

h : = [0 (r0~sin(S~) s) ( -ffl0r&~(l + w)sin(8~)s) 

- (r0~ + r0~ cos(8~) s) (©0 r0~ + ©0 r0~ ( 1 + m) cos(5~) e) 0] 


> ne:=simpiiTy(tayior(h[2], epsilon)); 

he := 

- ©0 r0~ 2 + 1 -2 r0~ 2 cos(8~) ©0 - r0~ 2 cos(8~) ©0 u\ e + 1-©0 r0~ 2 - r0~ 2 ©0 u\ e 2 




11 0 0 

r2 r0~ + r0~ cos(8~) s + - r0~ - - r0~ cos(8~) z + 


j j 3 5 

- - r0~ cos(8~) + - r0~ cos(S~) 3 s 3 + - r0~ cos(8~) 2 - l r0~ cos(8~) 4 - - r0~ s 4 + 

2 2 4 8 8 


5 7 3 

- r0~ cos(8~) 3 + - r0~ cos(8~) 5 + - r0~ cos(8~) e 5 + o|e 6 | 



e:=(evalm(-(r/r2Hcrossprod(h,v)/mu))): 


> ei :=simpnty(tayior(e[i ],epsilon,4)); 


3 

el := (3 cos(8~) + 2 cos(8~) u) e + - cos(8~) 2 + 3 u cos(8~) 2 + cos(8~) 2 u 2 + - + u s 2 

^ 2 

+ |cos(8~) 3 + 2 cos(8~) u + cos(8~) w 2 | e 3 + o|e 4 | 
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e2 := 0 


> e3:=simplity(taylor(e[3], epsilon, 4)); 


e3 := sin(8~) u z + 13 cos(5~) sin(5~) + 3 cos(8~) sin(5~) u + cos(5~) sin(5~) u 2 i e 2 + 

3 


3 3 

- sin( 8~ ) — - sin( 8~ ) cos( 8~ ) 2 + 2 sin( 5 ~)u + u 2 sin( 8~ ) 


8" +0|s 4 ' 


EXCENTRIOTY 

> ee:=simplify(taylor(sqrt(dotprod(e,e, 'orthogonal')), epsilon, 5)); 


ee 


:= J%Ts + ^cos(8~)j 


9 cos(S ~) 2 + 18 u cos(S ~) 2 + 9 + 18 u+ 12 cos( 8~) 2 w 2 + 2 w 3 cos(S ~) 2 + 10 u 2 + 2 u 3 
|/ J%1 z 2 -- ]%l\- 528 cos( 8~) 4 w 4 - 8 w 5 - 20 m 4 - 4 cos( 8~) 4 u 6 - 24 « 5 cos(8~) 2 

' O 

- 80 w 5 cos(S ~) 4 - 120 cos(S ~) 2 w 4 - 324 cos( 8~) 4 + 108 cos(S ~) 6 + 4 cos( 8~) 6 a 6 

- 2067 u 2 cos(S ~) 4 + 636 cos( 8 ~)^ u 2 + 48 cos( 8 ~)^ a 3 - 1552 a 2 cos( 8~) 4 
+ 252 cos( 8 ~)^ a 4 + 753 cos( 8 ~)^ u 2 + 432 cos( 8 ~)^ a - 1296 a cos( 8~) 4 

- 45 C 9 s( 8~) 2 u 2 - 24 a 3 - 9 u 2 - 132 a 2 cos(8~) 2 |/|81 cos(S ~) 4 + 216 a cos( 8~) 4 
+ 198 u 2 cos ( 8~) 4 + 18 cos( 8~) 2 u 2 + 72 a 3 cos( 8~) 4 + 24 u 3 cos( 8~) 2 

+ 9 cos( 8~) 4 « 4 + 6 cos( 8~) 2 a 4 + a 4 js 3 + o|e 4 j 

%1 := 9 cos( 8~) 2 +12 a cos( 8~) 2 + 3 cos( 8~) 2 u 2 + a 2 

> eel :=simplity(coert(ee, epsilon, 1 ),trig); 

eel := j 9 cos( 8~) 2 +12 a cos( 8~) 2 + 3 cos( 8~) 2 u 2 + u 2 

> ee 2 :=simplify(eval(coett(ee, epsilon, 2 ))); 

ee2 := - cos( 8 ~) | 


9 cos(8~) 2 + 18 a cos(8~) 2 + 9 + 18 a + 12 cos(8~) 2 u 2 + 2 u 3 cos(8~) 2 + 10 a 2 + 2 u 3 

| j J 9 cos(S~) 2 + 12 a cos(8~) 2 + 3 cos(8~) 2 u 2 + a 2 
ENERGY CONSTANT 

> EN:=simplify(taylor(dotprod(v,v,’orthogonar)/2-mu/r2, epsilon, 4)); 

2 + 2 a + a 2 | 2 
z + 

-3 + 5 cos( 8 ~) 2 j ^3 _ ^4 
2 r0~ 1 

SEMIMAJOR AXIS 


1 ^ cos(8~) 


r>r 1 l 1 , ^ cos(8~) (2 + a) 1 p. -3 cos(8~) + 

-2^: + 75: e+ 2 To: 
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> a:=simplify(taylor(-mu/EN/2, epsilon, 4)); 

a := r0~ + (4 r0~ cos(5~) + 2 r0~ cos(5~) u) e + 

1 13 r0~ cos(8~) 2 + 2 r0~ + 2 r0~ u + r0~ u 2 + 16 r0~ cos(8~) 2 w + 4 r0~ cos(8~) 2 u 2 \ 
2 I ^ 

s + ! 13 r0~ cos(8~) + 45 r0~ cos(b~y + 24 r0~ cos(8~) u + 16 r0~ cos(8~) « 2 
+ 84 r0~ cos(8~) 3 u + 4 r0~ cos(8~) u 3 + 48 r0~ cos(8~) 3 u 2 + 8 r0~ cos(8~) 3 u 3 e 3 
+ o|e 4 | 

> simplity(coert(a,epsiion,2)/r0); - — 


13 cos(8~) +2 + 2 u + u + 16 «cos(8~) + 4cos(8~) 2 w 2 

MEAIN ANGULAR RATE (TIMES CIRCULAR RATE) 

> n:=simplify(taylor(sqrt(rO*rO*rO/a/a/a), epsilon, 4)); 

n:= 1 + (-6 cos(8~) - 3 cos(8~) u) e + 


21 


cos(8~) 2 -3-3«--w 2 +6w cos(8~) 2 + ^cos(8~) 2 u 2 
Z Z 2 


2 o 

8 + 6cos(8~)m z 


+ 3 cos(S~) 3 w 2 + - cos(S~) 3 u + 9 cos(8~) u- ^cos(8~) 3 +- cos(8~) i ? 
•4 2 2 

3 


1 3 3 21 

+ - « cos(8~) + — cos(8~) 


8“+C^£ 4 | 


ORBITAL EEKIOU (TIMES 2Pi7CIRC'ULAR RATE) 
> t:=simplify(taylor((a/r0)**(3/2), epsilon, 3)); 


t 1 + (6 cos(8~) + 3 cos(8~) u) e + 

— cos(8~) 2 + 3 + 3 w + -« 2 + 30h cos(8~) 2 + -^cos(S~) 2 u 2 e 2 + o|e 3 | 


IMCKLASE IN PER1C EE~HEIGHT O VER CIRCULAR ORBI T 
> dhp:=simplify(taylor(a*(1-ee)-r0, epsilon, 3)); 


dhp :~\-rO~J%\ + 4r0~cos(8~) + 2r0~cos(8~)M|e-^r0~|81 cos(S~) 3 

+ 150 cos(8~) 3 u + 9 cos(8~)+ 18 cos(8~) u + 84 cos(8~) 3 u 2 + 14 m 3 cos( 8~) 3 
+ 18 cos(8~) u 2 + 6 cos(8~) w 3 - 26 cos(8~) 2 J%1 -4 J%1 - 4 u J%1 - 2 u 2 J%1 
-32cos(8~) 2 wy%T-8cos(8~) 2 u 2 J%\ \j J%1 e 2 + o|e 3 | 

%1 9 cos(8~) 2 +12 u cos(8~) 2 + 3 cos(8~) 2 u 2 + u 2 

INCREASE IN APOCEE HEIGHT 

> dha:=simplify(taylor(a*(1+ee)-r0, epsilon, 3)); 


dha r0~ J% 1 + 4 r0~ cos(8~) + 2 r0~ cos(8~) u \ 8 + — — 1 8 1 cos(8~) 3 

+ 150 cos(8~) 3 u + 9 cos(8~) + 18 cos(8~) u + 84 cos(8~) 3 u 2 + 14 w 3 cos(S~) 3 
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+ 18 cos(5~) u 2 + 6 cos(5~) h 3 + 26 cos(5~) 2 J%T + 4 J%\ +4 u J%1 +2u 2 j%\ 
+ 32 cos(5~) 2 u J%\ + 8 cos(5~) 2 u 2 j%l\/ ]%l s 2 -f- o|e 3 | 

%1 := 9 cos(6~) 2 + 12 u cos(5~) 2 + 3 cos(8~) 2 u 2 + u 2 

> dhpt:=convert(simplity(subs({u=sqrt(3*sin(theta’Fi/l8)”2-3 w sin(delta)’*2),ru=biS/a.i6} 1 dhp)),polynom 
>): 

> dhpts:=subs(epsilon=U.9l46’20/6378.l6,dript) $ theta=3..6: 

DELI A INFLUENCE IN PERIGEE HEIGHT FOR !HETAM=J0-60 DEG (KM) 

> plot({dhpts},delta=-Pi/6..Pi/6); 



> dhat:=convert(simplity(subs({u=sqrt(3’sin(theta’Pi/l8)’"2-3 w sin(delta)”2),rt)=637a.i6},dha)) 1 polynom 
>): 

> dhats:=subs(epsilon=U.9l4b 


Appendix E: DECAY1.MS 


INITIAL TRUE ANOMALY 

> sinnuO:=simplify(taylor(dotprod([ 0 , 1 ,0],crossprod(r 1 e),'orthogonar)/r2/ee, epsilon, 5 )); 
sin(5~)w 1 . 

sumuO \= t= sm(5~) u cos(5~) 

J% 1 2 ' 

|9 cos(8~) 2 + 24 u cos(8~) 2 + 18 cos(8~) 2 u 2 -8u 2 + 4u 3 cos(5~) 2 -9-lSu / 

3/2 1 I * - 

%1 s + - sin(6~) w|378 cos(8~)° - 792 u 3 cos(5~) 2 + 1836 cos(S~) 6 u 

+ 3171 cos(5~) 6 w 2 - 999 cos(8~) 2 u 2 + 24 u 3 + 16 w 4 + 9 u 2 + 176 u 5 cos(S~) 4 

- 324 cos(5~) 2 w 4 + 12 cos(S~)^ - 21 « 2 cos(8~) 4 + 468 cos(8~) 4 u 4 
+ 20 cos(S~) 4 u 6 - 48 w 5 cos(5~) 2 + 1032 cos(8~) 6 u 4 - 108 u cos(8~) 4 

- 162 cos(5~) 2 - 648 u cos(8~) 2 + 2568 cos(6~) 6 w 3 + 192 cos(5~) 6 u 5 

+ 400 w 3 cos(6~) 4 J / %1 e 2 + o(g 3 | 


%1 9 cos(8~) 2 + 12 u cos(5~) 2 + 3 cos(8~) 2 m 2 + 1 / 2 

nu(J:=simpnty(simpiity(tayior(arcsin(sinnuU), epsilon, b),sqrt,assu me=positive)); 


vO := - arcsin 


sin(8~) u 


j%i 


- 2 sin(5~) u 


|9 cos(8~) 2 + 24 u cos(8~) 2 + 18 cos(8~) 2 u 2 -8u 2 + 4« 3 cos(8~) 2 - 9 - 18 u\ 

a 


signum(cos(8~)) 


%1 


3/2 |(3 + 2 w) 2 


%1 


£ + — j 8 cos(8~) 4 + 8 cos(8~)' 


+ 68 u 5 cos(S~) 2 - 16 « 5 + 108 m 5 cos(8~) 4 - 124 u 4 + 180 cos(S~) 2 u 4 
+ 540 cos(8~) 4 u 4 - 354 m 3 + 152 w 3 cos(8~) 2 + 1302 u 3 cos(8~) 4 -486 w 2 
- 24 cos(S~) 2 u 2 + 1590 w 2 cos(8~) 4 + 918 1 / cos(8~) 4 - 324 u - 54 u cos(8~) 2 - 81 


+ 189 cos(8~) 4 jsm(8~) « signum(cos(8~)) cos(8~)/ %1 

%1 := 9 cos(8~) 2 + 12 m cos(8~) 2 + 3 cos(8~) 2 w 2 + m 2 
> simplify(coett(nuO, epsilon, 0)); 


5/2 j(3+2 u) 2 


%1 


s 2 + 0 


-arcsin 

sin( 8~)u 


j 9 cos(8~) 2 +12 u cos(8~) 2 + 3 cos(8~) 2 u 2 + « 2 


^ simpniy isimpinyicoemnuu, epsilon, 1 ),sqrt,assume-positive),pow); 
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- ^ signum(cos(8~)) |2 cos(8~) 2 «^-4« + 6u cos(8~) 2 -3 + 3 cos(8~) 2 j u sin(8~) 


1 


j 9 cos(8~) 2 +12 u cos(8~) 2 + 3 cos(8~) 2 u 2 + w 2 


|9 cos(8~) 2 + 12 u cos(S~) 2 + 3 cos(8~) 2 u 2 + « 2 


3/2 


> factor(4~u A 3+i8 w u A 2+24*u+9); 

(3 + 2w)|2w 2 + 6w + 3| 

> factor{8*u A 2+1 B*u+y); 

INITIAL ANGULAR SPEED OF 1 

(3 + 2 m) (4 « + 3) 

ELLIPTIC ORBITAL FRAME — 


NOTE MAPLEV'S GLORIOUS STUPIDITY IN EASY SIMPLIFICATIONS, 

WHILE HANDLING DIFFICULT ONES WELL. 

> nup:=subs(signum(cos(delta))=1,simplify(taylor(n*(1+ee*cos(nuO)) A 2/(1-ee*ee) A (3/2),epsilon t 3))); 


nup 1 +i2 cos(8~) 


(3 + 2ay 


J 


%1 


- 6 cos(8~) — 3 cos(8~) u 


s - 


-270 u cos(8~y 


- 216 cos(8~) 2 u 2 - 9 u 2 + 82 cos(8~) 2 w 4 - 81 cos(8~) 2 + 20 u 5 cos(8~) 2 
+ 60 m 3 cos(8~) 4 - 18 w 3 + 30 w 3 cos(8~) 2 + 558 cos(8~) 4 w 4 - 8 w 4 
+ 1992 u 2 cos(8~) 4 + 341 1 u 2 cos(8~) 4 + 2808 u cos(5~) 4 + 891 cos(8~) 4 


+ 3 u !%\ 


(3 + 2 “ r +3 l%[ K 3 + 2w > -2 


%\ 


%1 


« -297cos(8~r /%1 


(3 + 2«V 


J 


%1 


+ 27 cos( 8~) 2 J%T \ i3+ ^ u) --30u 4 J%l 

J 


< 3 + 2 "> co S <8~) 4 


%1 


- 10 « 4 7%1 j' 3 ^ — cos (8~) 2 -720 uj%\ cos(8~) 4 

+ 63 uj%l j ° cos(8~) 2 + 12 u 2 J%[ j (3+ ^, U)2 cos(8~) 2 


- 27 a 3 J%T j ( - cos(5-) 2 -62\u 2 J%1 cos(8-) 4 


%1 


228 k 3 J%1 — + J U ^ cos(8~) 4 

g %1 


%1 


3/2 (3 + 2u) 2 


%1 


s 2 + C^s 31 
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%1 := 9 cos(5~) + 12 u cos(8~) 2 + 3 cos(8~) 2 u 2 + w 2 

WITH A HTTEE COAXING, IT WILL DOIT: 

> nupl :=simplify(subs(9*cos(delta) A 2+12*cos(delta) A 2‘u+3*cos(delta) A 2*u A 2+u A 2=c,coeff(nup epsilon 

> 1)l.sart.assume=DositiveV r 

nupl := cos(5~) « 


1)),sqrt,assume=positive); 

> nup2:=coeff(nup, epsilon, 2); 


nup2 := - 


-270 u cos(5~) 2 - 216 cos(8~) 2 « 2 - 9 u 2 + 82 cos(8~) 2 « 4 - 81 cos(8~) 2 

+ 20 w 5 cos(8~) 2 + 60 u 5 cos(8~) 4 - 18 u 3 + 30 u 3 cos(8~) 2 + 558 cos(8~) 4 « 4 - 8 « 4 
+ 1992 w 3 cos(8~) 4 + 3411 u 2 cos( 8~ ) 4 + 2808 wcos(8~) 4 + 891 cos(8~) 4 

+ 27 cos(8~) 2 J%T j ^ U)2 - 30 m 4 J% I 


(3+2 "> ,4 

cos(S~) 


J %1 


- 10 a 4 J%1 j (3+ 0/ ^ U) cos(8~) 2 - 720 cos(8~) 4 


63 u J% 1 


'(3 + 2uy 
%1 


cos(8~) 2 + 12 m 2 J%1 ■ — - + - — /-s \2 


J %1 


cos(8~)' 


- 27 w 3 y%T j cos(8~) 2 - 621 w 2 J%T j cos(8~) 4 


- 228 w 3 J % T 


(3 + 2 "> cos (5 -) 4 


%1 


%1 


3/2 ; (3 + 2 m) 2 


%1 


%1 := 9 cos(8~) 2 + 12 a cos(8~) 2 + 3 cos(8~) 2 a 2 + a 2 

> nup^s:=suPs({(a'cos(aeitaj“z-*-i^~cos(delta) A 2*u+3*c os(delta) A 2*u A 2-i-u A 2)=c,sqrt((3+2*u) A 2/(9*cos( d 

> elta) A 2+12*cos(delta) A 2*u+3*cos(delta) A 2*u A 2+u A 2))=(3+2*u)/sqrt(c)},nup2); 

:= - [-270 u cos(8~) 2 - 216 cos(8~) 2 u 2 - 9 a 2 - 621 u 2 (3 + 2 u) cos(S~) 4 

- 228 a 3 (3 + 2 a) cos(8~) 4 - 27 a 3 (3 + 2 a) cos(8~) 2 - 30 a 4 (3 + 2 a) cos(8~) 4 

- 10 a 4 (3 + 2 a) cos(8~) 2 + 27 cos(8~) 2 (3 + 2 u) - 720 a (3 + 2 u) cos(8~) 4 

+ 63 u (3 + 2 a) cos(8~) 2 + 12 u 2 (3 + 2 u) cos(8~) 2 + 3 a 3 (3 + 2 u) + 3 (3 +2 u) u 2 

- 297 cos(8~) 4 (3 + 2 a) + 82 cos(8~) 2 a 4 - 81 cos(8~) 2 + 20 a 5 cos(8~) 2 
+ 60 a 5 cos(8~) 4 - 18 a 3 + 30 a 3 cos(8~) 2 + 558 cos(8~) 4 a 4 - 8 a 4 
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3/2 

+ 1992 u r cos(8~)^ + 341 1 u ^ cos(5~)^ + 2808 u cos(5~)^ + 891 cos(S~)^j/jc 
j ( 3+2«) 2 

J c 

> nup2ss:=simplify(simplity(nup2s,sqrt) ,, (c A (372)’sqrt((3+2’u) / '2/c))/c/(3+2*u), radical); 

nup2ss := - j-cos(8~)^ P + 6 P cos(8~)^ - P - 12 u cos(8~)^ + 24 u cos(8~)^ 

+ 18 cos( 8~)4 - 9 cos(8~)^| u/c 

> mjp2sss:=subs(c=9*cos(delfa) A 2+1Z’cos(delta) A 2’u+3’cos(delta) A 2’u A 2+u A 2,nup2ss); 

nup2sss := - |-cos(8~)^ P + 6 P cos(8~)^ - P - 12 u cos(8~)^ + 24 u cos(8~)^ 

+ 18 cos( 8~)4 - 9 cos(8~)^| u j\9 cos(8~)^ + 12 u cos(8~)^ + 3 cos(8~)^ u ^ + i P 

FINALLY: ~ 

> nup22:=simplify(nup2sss, radical); 


nup22 := -u 2 cos(8~) - 1 

WE’LL CHECK FOR SIMPLIFICATION ERRORS: O.K. 

> plot(subs(cos(delta)=0.5,nup2-nup2sss),u=-1.5..1.5); 


> nudiffO:=l +nup'Tepsiion+nup22’ r epsiion A 2; 


nudiffO := 1 + cos(S~) ue-u 2 cos(8~)^ -1 ^ 
> psiO:=simplity(taylor((u+l -nuditt0)/nudiftu,epsilon,3)); 


| 21 

y/0 := u + i-cos(8~) u - cos(8~) u [e + 

|2 u cos(8~)^ - u + 3 cos(8~)^ P -p + u 3 
> tactor(coett(psiU, epsilon, 2)); 


cos(8~)‘‘ 


s 2 + OI 


m(1 + m) «cos(8~)‘ 3 + 2cos(5~)^- 1 

ANGULAR SPEED IN C.O.M: ELLIPTIC ORBIT 

> nudiff:=subs(signum(cos(delta))=1,simplify(taylor(n*(1+ee*cos(nu)) A 2/(1-ee*ee) A (3/2), epsilon, 4))); 

nudiff'.- 1 + J 2 J% 1 cos(v) - 6 cos(8~) - 3 cos(8~) w|e — 13 u J% 1 

+ 99 cos(8~) 3 cos(v) - 24 cos(8~)^ J% 1 + 3 J% 1 -24 cos(8~)^ w J% 1 

- 6 cos(8~)^ J%1 - cos(v)^ 7%T - 12 cos(v)^ J%1 u cos(8~)^ 

- 3 cos(v)^ J% 1 cos(8~)^ P — 9 cos(v)^ J%1 cos(8~)^ + 2 cos(8~) cos(v) 

3 3 2 

- 18 cos(8~) cos(v) u + 4 cos(8~) cos(v) u + 96 cos(8~) cos(v) u 

+ 16 cos(8~) 3 cos(v) w 3 + 180 cos(8~) 3 cos(v) u- 9 cos(8~) cos(v)[/ J% 1 P + ^ 


72 cos(8~) cos(v)^ 7%T m 3 + 36 cos(8~) cos(v)^ J% 1 P 
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+ 1080 cos(8~) 3 cos(v) 2 J%1 u - 3564 cos(v) cos(8~) 4 - 1416 cos(v) u 3 cos(8~) 2 

- 13284 cos(v) u cos(8~) 4 - 603 cos(v) cos(8~) 2 u 2 + 864 cos(8~) 3 J% 1 

+ 9 cos(v) u 2 + 24 cos(v) u 3 - 16 cos(v) u 5 - 4 cos(v) u 4 + 144 cos(8~) u 3 J%\ 

+ 96 cos(8~) u 2 J%1 + 2448 cos(8~) 3 u J%1 + 324 cos(8~) 3 cos(v) 2 J%1 
+ 52176 cos(8~) 6 cos(v) u 3 + 1104 cos(8~) 4 cos(v) u 4 + 20772 cos(S~) 6 cos(v) u 4 
+ 24 cos(8~) 2 cos(v) + 356 cos(8~)^ cos(v) u ^ + 196 cos(8~) 4 cos(v) u ^ 

+ 4272 cos(8~) 6 cos(v) u 5 + 1184 cos(8~) 4 cos(v) « 5 + 49140 cos(8~) 6 cos(v) u 

- 120 cos(8~) 5 J% 1 w 5 - 840 cos(v) 2 cos(8~) 5 J%1 u 4 - 8880 cos(8~) 5 J%\ w 2 

- 96 cos(8~) 3 y%T « 4 + 48 cos(8~) J%1 u 4 + 70959 cos(8~) 6 cos(v) u 2 
+ 16 cos(v) 2 cos(8~) y%T « 4 - 4 cos( v) 2 cos(8~) J%1 u 5 

+ 48 cos(v) 2 cos(8~) 3 J%1 u 3 - 40 cos(v) 2 cos(8~) 3 J%1 u 5 

- 84 cos(v) 2 cos(8~) 5 J%1 u 5 - 5940 cos(v) 2 cos(8~) 5 J%1 u 2 

- 168 cos(v) 2 cos(8~) 3 J%1 u 4 + 13500 cos(8~) 6 cos(v) 

- 3240 cos(v) 2 cos(8~) 5 J%1 u 3 - 5076 cos(v) 2 cos(8~) 5 j"%\ u 

- 4680 cos(8~) 5 J%\ u 3 - 8160 cos(8~) 5 J%T u + 2128 cos(8~) 3 J%1 u 2 

- 40 cos(5~) 3 J%\ u 5 + 936 cos(v) 2 cos(8~) 3 J%1 u 2 - 15861 cos(8~) 4 cos(v) u 2 

- 6008 cos(8~) 4 cos(v) u 3 - 96 cos(S~) 2 cos(v) u 5 - 864 cos(8~) 2 cos(v) u 4 

- 1620 cos(v) 2 cos(8~) 3 y%T-2880 cos(8~) 3 J%1 +528 cos(8~) 3 J%T u 3 

3/2 

-I200cos(8~) 5 y%r u 4 |/%1 8 3 + o|e 4 ! 

%1 := 9 cos(S~) 2 +12 u cos(8~) 2 + 3 cos(8~) 2 u 2 + u 2 

SlINCK €OEFFTCENTS SEEM lOTJROW, WE WILL CHECK TH AT CHEATER 

ORDER TERMS REMAIN SMALL, GIVING VALUES TO THETA AND DELTA 
FOR A COMPLETE ORBIT IN NU: 

> nudiffs:=convert(nudiff,polynom); 

nudiffs := 1 +[2^%! cos(v) - 6 cos(8~) - 3 cos(8~) u\z -[3 u J%\ 

+ 99 cos(8~) 3 cos(v) - 24 cos(8~) 2 J%\ + 3 J%1 - 24 cos(8~) 2 u J%T 

- 6 cos(8~) 2 u 2 J%1 - cos(v) 2 y%T u 2 - 12 cos(v) 2 J%1 u cos(8~) 2 

- 3 cos(v) 2 J% 1 cos(8~) 2 u 2 — 9 cos(v) 2 J% 1 cos(8~) 2 + 2 cos(8~) cos(v) u 2 

- 18 cos(8~) cos(v) u + 4 cos(8~) cos(v) u 3 + 96 cos(8~) 3 cos(v) u 2 
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+ 16 cos(8~) 2 cos(v) u 3 + 180 cos(8~) 2 cos(v) u - 9 cos(8~) cos(v)| e 2 / J % I + -^ j 
72 cos(8~) cos(v) 2 /% 1 m 2 + 36 cos(8~) cos(v) 2 J% 1 h 2 

+ 1080 cos(8~) 2 cos(v) 2 J%1 u - 3564 cos(v) cos(8~) 4 - 1416 cos(v) w 3 cos(8~) 2 

- 13284 cos(v) u cos(8~) 4 - 603 cos(v) cos(8~) 2 u 2 + 864 cos(8~) 2 J% 1 

+ 9 cos(v) u 2 + 24 cos(v) u 3 - 16 cos(v) m 2 - 4 cos(v) u 4 + 144 cos(5~) u 3 J%1 
+ 96 cos(5~) u 2 j%\ + 2448 cos(8~) 2 u J%1 + 324 cos(8~) 3 cos(v) 2 J%1 
+ 52176 cos(5~)^ cos(v) u 3 + 1104 cos(8~) 4 cos(v) m 4 + 20772 cos(8~)^ cos(v) u 4 
+ 24 cos(8~) 2 cos(v) u ^ + 356 cos(8~)^ cos(v) u ^ + 196 cos(8~) 4 cos(v) it* 3 
+ 4272 cos(5~)^ cos(v) u 3 + 1184 cos(8~) 4 cos(v) u 3 +49140 cos(8~)^ cos(v) u 

- 120 cos(5~)^ /%1 u 3 - 840 cos(v) 2 cos(8~)^ J% 1 u 4 - 8880 cos(5~)^ J% 1 u 2 

- 96 cos(5~) 3 J%T m 4 + 48 cos(S~) J%\ u 4 + 70959 cos(5~) 6 cos(v) « 2 
+ 16 cos(v) 2 cos(5~) j% 1 m 4 - 4 cos(v) 2 cos(8~) J%1 

+ 48 cos(v) 2 cos(5~) 2 J% 1 w 2 - 40 cos(v) 2 cos(8~) 2 J% 1 m 2 

- 84 cos(v) 2 cos(S~)^ J%i u 3 - 5940 cos(v) 2 cos(8~)^ J%1 u 2 

- 168 cos(v) 2 cos(8~) 2 J%\ w 4 + 13500 cos(8~) 6 cos(v) 

- 3240 cos(v) 2 cos(8~) 2 j%\ u 3 - 5076 cos(v) 2 cos(8~)^ J%1 u 

- 4680 cos(8~) 5 J%\ u 3 - 8160 cos(8~) 5 J%\ u + 2128 cos(8~) 3 j%T u 2 

- 40 cos(S~) 3 J%\ u 5 + 936 cos(v) 2 cos(8~) 3 J%\ u 2 - 15861 cos(8~) 4 cos(v) u 2 

- 6008 cos(8~) 4 cos(v) u 3 - 96 cos(8~) 2 cos(v) u 3 - 864 cos(S~) 2 cos(v) w 4 

- 1620 cos(v) 2 cos(8~) 5 J%\ - 2880 cos(8~) 5 J%\ + 528 cos(8~) 2 J%\ u 3 

2/2 

- 1200 cos(8~) 5 J%1 u 4 \e 3 /%l 

%\ := 9 cos(S~) 2 + 12 u cos(8~) 2 + 3 cos(8~) 2 u 2 + u 2 
TERM OK ORDER EPSILON (DELTS“0;THETA nkADlAN): 

> nud1:=eval(subs({theta=1,delta=0},subs(u=sqrt(3*sin(theta) A 2-3*sin(delta) A 2),coeff(nudiffs, epsilon, 1) 

> ))); 

nudl := 2 J 9 + 12 JJ Jsin(l) 2 + 12 sin( 1 ) 2 cos( v) - 6 - 3 p j sin( 1 ) 2 
EPSILON SQUARE: 

> nud2:=eval(subs({theta=1,delta=0},subs(u=sqrt(3*sin(theta) A 2-3*sin(delta) A 2),coeff(nudiffs, epsilon, 2) 

> ))); 
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nud2 := - 


-21 T 2 J sin( 1 ) 2 y%T + 90 cos(v) - 21 J%\ - 18 sin( 1 ) 2 J%T 


- 12 cos( v) 2 7%T sin( 1 ) 2 - 12 cos(v) 2 J%\ J 3 J sin( 1 ) 2 - 9 cos(v) 2 f %1 


3/21 


+ 294 cos(v) sin( 1 ) 2 + 162 cos(v) JJ J sin( 1 ) 2 + 60 cos(v) JJ |sin( 1 ) 2 | j J%\ 
%1 := 9 + 12 J3 Jsin( 1 ) 2 + 12 sin( 1 ) 2 

~EPSIEDNXTJBE: 

> nud3:=eval(subs({theta=1,delta=0},subs(u=sqrt(3*sin(theta) A 2-3*sin(delta) A 2),coeff(nudiffs, epsilon, 3) 
^ )))» 


nud3 := 


-19968 sin(l) J% 1 + 163512 cos(v)sin(l) 2 


3/2 5/2 

-9360cos(v) 2 y%ry3 |sin(l) 2 | - 1152 cos(v) 2 y%F y3 |sin(l) 2 | -2016 /%T 


+ 9936 cos( v ) - 3996 cos( v ) 2 J3 J sin( 1 ) 2 + 35856 cos( v ) J 3 j sin( 1 ) 2 

- 1296 cos(v) 2 7%T - 5712 J3 Jsin(l) 2 j%T- 14904 cos(v) 2 7%T sin( 1 ) 2 

3/2 5/2 

+ 134328 cos(v) j3 |sin(l) 2 | + 48096 cos(v) yj|sin(l) 2 | 

3/2 

- 8928 cos(v) 2 J%\ sin( 1 ) 4 - 12024 Js |sin( 1 ) 2 | J%1 


5/2 


- 1440 J%T ^ sin( 1 ) 2 j + 189072 cos(v) sin( 1 ) 4 + 15552 cos(v) sin( 1 ) 6 


-11232 /% 1 sin(l) 


%1 


3/2 


%1 := 9 + 12 JJ Jsin( 1 ) 2 + 12 sin( 1 ) 2 

EES1LON BY ANOIHERTNAME, SO 1 HA 1' IT WI LL NOT BBT UKB I HE KEST OF 
THE CALCULATIONS: 

> sigma:=0.91 45*20/6378. 16; 

ct := .002867598179 

THE THIRD TERM GETS MUCH BIGGER THAN THE SECOND: 

> plot({nud2,nud3},nu=0..2*Pi); 



mgTgjTjj 


T 




* *rr’' 



•2*Pi); 



IE FIRST: 
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BU I EPSILON PUTS THINGS IN PERS PECTIVE: WE NEE1) ONLY KEEP TERMS 
TO ORDER EPSILON SQUARE. 

> plot({nud1 ,sigma*nud2},nu=0..2*Pi); 



AND NOW THE WHOLE PICT URE: 

> plot({1+sigma*nud1,sigma*sigma*nud2},nu=0..2*Pi); 
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WE WILL CONTINUE SIMPLIFYI NG THE TER MS OF INSTANT ANGU LAR SPEED 

> nudiff1:=coeff(nudiff, epsilon, 1); 

nudiffl := 

2 J9 cos(8~) 2 +12 u cos(8~) 2 + 3 cos(8~) 2 u 2 + u 2 cos(v) - 6 cos(8~) - 3 cos(8~) u 

> nudin 2 :=coert(nuditt,epsiion, 2 ); 

nudiff2 - 1 3 w ]%l + 99 cos(8~) 3 cos(v) - 24 cos(S~) 2 J%\ + 3 J%1 
-24 cos(8~) 2 u J%\ - 6 cos(8~) 2 u 2 J%\ - cos(v) 2 J%1 u 2 

- 12 cos(v) 2 J%1 u cos(8~) 2 - 3 cos(v) 2 ]%l cos(8~) 2 u 2 

- 9 cos(v) 2 /% 1 cos(8~) 2 + 2 cos(8~) cos(v) u 2 - 18 cos(8~) cos(v) u 

3 3 2 3 3 

+ 4 cos(S~) cos(v) u + 96 cos(8~) cos(v) u + 16 cos(8~) cos(v) w 

+ 180 cos(8~) 3 cos(v) u - 9 cos(S~) cos(v)|/ J%1 

%1 := 9 cos(8~) 2 +12 u cos(8~) 2 + 3 cos(8~) 2 u 2 + u 2 

> nuditT 2 :=collect(simplify(nudin 2 ,radical,assume=positive) ) cos(nu)); 


nudiff2 := 


-9 cos(8~) 2 y%r-M 2 y%T-12cos(8~) 2 M 7%T-3cos(8~) 2 m 2 J%T!cos(v) 2 


J%\ 

99 cos(8~) 3 + 2 cos(8~) w 2 - 18 cos(8~) m + 4 cos(8~) w 3 + 96 cos(8~) 3 u 2 
+ 16 m 3 cos(8~) 3 + 180 cos(8~) 3 u-9 cos(8~)|cos(v)/ J% 1 
3 « y%T - 6 cos(8~) 2 m 2 7%T - 24 cos(8~) 2 J%T + 3 J%T - 24 cos(8~) 2 u J%T 

7 ^ 
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%1 := 9 cos(8~) 2 + 12 u cos(S~) 2 + 3 cos(8~) 2 w 2 + m 2 

> cO:=simplify(coett(nuditr2,cos(nu),U), radical, assume=positive); 

cO := -3 u + 6 cos(8~) 2 u 2 + 24 cos(8~) 2 -3 + 24 u cos(8~) 2 

> factor(coeff(c0,cos(delta)**2)); - — 

6 (2 + w ) 2 

> factor(coeff(c(J,cos(delta),0)); — — 

-3 m - 3 

> cl:=collect(coeft(nudift2,cos(nu),l),cos(delta)); 

|99 + 96 w 2 + 16 w 3 + 180 i/|cos(8~) 3 

J 9 cos(8~) 2 + 12 u cos(S~) 2 + 3 cos(S~) 2 u 2 + w 2 
[2 m 2 - 18m + 4m 3 -9|cos(S~) 

J 9 cos(8~) 2 + 12 u cos(S~) 2 + 3 cos(5~) 2 w 2 + u 2 

> tactor( 1 6*u A 3+yy+ 1 8U*u+96*u A 2); 

> factor(-9+4*u A 3-1 8*u+2*u A 2); 

99 + 96 h 2 + 16m 3 + 180 u 
(2 m + 1 ) 1 2 w 2 — 9| 

> c2:=simplify(coeff(nucJitr2,cos(nu),2),assume=positive); “ 

c2 := 9 cos(8~) 2 + 12 u cos(8~) 2 + 3 cos(8~) 2 w 2 + « 2 

> factor(coett(c2,cos(delta)”*2)); 

3 (w + 3)(l +u) 
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DECAY2.MS 

DECAY MINIMIZATION: FINAL SATELLITE ORBIT 

> r:=array([r0'(1+rlx'epsHon),0,r0’epsilon’rlz]); 

r := [rO ( 1 + rix s) 0 rOzrlz ] 

> v:=array([omega0*r0*epsilon*v1x,0,omega0*r0*(1+v1z*epsilon)]); 

v := [©0 rO 8 vlx 0 ©Or0(l + vizs)] 


> with(linalg): 

> h:=crossprod(r,v); 


A:=|0 r0 1 2 e 2 rlz toO vlx — rfP" ( 1 + rix e) ©0 ( 1 + viz e) 0 


> h[2]:=collect(simplify(h[2]), epsilon); 

h ^ := | r0 2 rlz ©0 vlx - rO 3 ©0 rlx vlz\ e 2 + j -r0 2 ©0 viz - r(P“ ©0 rlx J s - rO ^ ©0 

> h2:=taylor(simplify(eval(h[2]*h[2])), epsilon); — 

h2 := rO ^ © 0 2 + r0 ‘ ^ © 0 2 (2 viz + 2 rlx) s + 
r0‘ ^ ©0 2 1 -2 rlz vlx + 2 rlx viz + (viz + rix) 2 | s 2 + 

2 rO ^ ©0 2 ( viz + rix) (-rlz vlx + rlx viz ) s 2 + rO ^ ©0 2 (-rlz vlx + rlx viz 


> resc:=simpiity(taylor(sqrt(dotprod(r,r, 'orthogonal’)), epsilon, 4), sqrt, assume=positive); 

resc := rO + rO rlx z + \-r0 rlz 2 e 2 - ^ rO rlx rlz 2 s 2 + ol I 

2 2 11 

> omega0:=sqrt(mu/r0/r0/r0); 


S 

0 


\ 

^ — ■— art : nr m — m — _i — z -1/ t _ in iw — 

rO 3 


> EN:=simplify(taylor(dotprod(v,v,’orthogonal')/Z-mu/resc, epsilon, 4), sqrt); 

1 u \x(vlz + rlx) 1 u viz 2 + vlx 3, + rlz^ - 2 rix 2 j 7 

2 rO rO 2 rO 

1 n rlx | -3 riz 2 + 2 rix 2 | 3 J 4 

2 r0 £ +C ^ £ 

> as:=simplify(tayior(-mu/EN/2, epsilon, 4)); 

ay := rO + 2 (viz + rix) r0 8 + (5 viz 2 + vix 2 + riz 2 + 2 rix 2 + 8 rlx vlz\ rO s 2 + |rix riz 2 


+ 2 rix 2 + 12 viz 2 + 4 viz vix 2 + 4 viz riz 2 +16 viz rix 2 + 28 rix viz^ + 4 rix vix 2 j 

r0 s 2 + ols^l 


> evs:=evalm(-r/resc-crossprod(h,v)/mu); 

r0( 1 +rix e) 


ev5 := 


1 0 ? 1 9 3 14 

rO + rO rlx e + -r0 rlz A e - - rO rlx rlz z e + Os 

2 2 1 
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rO 2 rlz 


rO- 


vlx - rO * 


| , rlx viz 

s 2 + 

-rtf 2 j 

rO 3 


J 


rO' 


vJz - rO 


2 H 


rO- 


rlx 


-rO 


2 H 


rO- 


rO~ 


rO{\ + v/zs)/p. 0 


rO e rlz 


rO + rOrlx 8 + ^ rO rlz 2 e 2 - ^ rO rlx rlz 2 e 3 + o|e 4 | 


rO 2 rlz 


rO- 


vlx - rO 1 


1 r/x v/z 

c 2 + 

| -rtf 2 I 

j rO 3 


i 


-rO 2 JL Wz _ r0 2 JL rIx 


rO- 


rO~ 


- ri r 


K 

| - rOevlx/n 

JrO 3 

JrO 3 


> es:=simplify(taylor(sqrnaotproa(evs,evs, 'orthogonal')), epsilon, 4), sqrt); 

es := J4 viz 2 + 4 rlx viz + r/x 2 + r/z 2 + 2 r/z v/x + v/x 2 e + ^(4 v/z 3 + 10 rlx viz 2 

+ 2 viz rlz 2 - 2 viz rlz vlx + 4 viz rlx 2 - 2 vlx rlx rlz - rlx rlz 2 + 2 viz vlx 2 

+ 2 rlx vlx 2 1 j J 4 viz 2 + 4 r/x viz + rlx 2 + rlz 2 + 2 rlz vlx + vlx 2 e 2 + o|e 3 
ANOTHER WAY TOX'OMPUTITEXL'LN TRIC'ITY: — 


es2:=simplify(1+2*h2*EN/mu/mu); 


es2 := 


2 \r0 m + 2 rO n (viz + rlx) e + rO n |-2 r/z v/x + 4 r/x v/z + viz 2 + r/x 2 | s 2 + 
2 rO ji ( viz + rlx) (-rlz vlx + rlx viz) e 3 + rO p (-rlz vlx + rlx viz) 2 s 4 


-IJU 

2 rO 


(v/z + r/x) v/z 2 + v/x 2 + r/z 2 - 2 r/x 2 2 1 H rlx -3 r/z 2 + 2 r/x 2 1 7 

£ + ~ — £ -I L o J 



rO 

&t 2 

rO 

C + 2 

rO 

8 + 

J 41 

± 2 

/ 2 
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I* 






> esz:=simpiify(taylor(esz, epsilon)); 

es2 := 1 4 v/z 2 + 4 r/x v/z + r/x 2 + r/z 2 + 2 r/z v/x + v/x 2 | e 2 + 14 v/z 3 + 10 r/x v/z 2 
+ 2 v/z r/z 2 - 2 v/z r/z v/x + 4 v/z r/x 2 - 2 v/x r/x r/z - r/x r/z 2 + 2 v/z v/x 2 
+ 2r/xv/x 2 je 3 + o|e 4 
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DECAY3.MS 

ORBITAL DECAY RATE: DECREASE IN SEMIAXIS PER REVOLUTION. 

> mu:=3.98618e5; 

> cosi:=evalf(cos(57*Pi/180)); 

> h:=59.77; 

> we:=7.292115e-5; 

H := 398618. 
cosi := .5446390348 
h := 59.77 


we:= .00007292115 

> K:=evalt(2*Pi/1 1 8*1 e-6*3^2); 


:= .1757161993 10‘ 8 

> p:=K’a A Z’'exp((6678.l6-a)/h)"(l-2’we’sqrt(a'a"a/mu)’cosi)’(Besseli(0,a'e/h)+2'e*Bessell(l,a*e/h 

> )); 


, 8 7 (11 1.7309687 -.01673080141 a) 

p := .1757161993 10 6 a z e 1 1 - .1258096100 10 

(Bessell(0, .01673080141 ae) + 2e Bessell( 1, .01673080141 a e )) 


-6 / 3 


a 


> with(plots): ~ ~ ~ 

> setoptions3d(axes=NORMAL,style=PATCHCONTOUR,contours=12,shading=Z): 

> plot3d(p,a=6675..6878,e=0..0.03,title='Semiaxis Decay'); 

Semiaxis Decay 



THE MINIMUM DECAY OCCURSTOTTMAXIMLIM a ANLTM1N1MUM e: 
> pe:=subs(e=ee/300,evalf(p)); 
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.1757161993 ICT 8 e ( ' ' '- 6807763 ' 01673080141 a) 


1. -.1258096100 10 
+ .006666666666 ee 


J ' (Bessell(0, .00005576933803 a ee 
Bessell( 1, .00005576933803 a ee)) 


) 


> pde:=evalf(pe) $ ee=0..6: 

> plot({pde),a=6678..7000,title= Decay Variation with Excentricity: e=0-0.02'); 


Decay Variation with Excentricity: e=0-0.02 



DECAY DECREASES WITH CKOWING a AMU WITH DIMINISHING e 

> pa:=subs(a=6378+50*i,p): 

> pda:=pa $ i=6..10: 

> plot({pda},e=0.. 0.02, title=' Decay variation with excentricity, a=300-500 Km ); 


G 2 
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Decay variation with excentricity, a=300-500 Km 
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DECAY4.MS 

DECAY RATE SIMPLIFIED EXPRESSION 

> mu:=3.9861 8e5; 

> cosi:=evalf(cos(57*Pi/1 80)); 

> h:=59.77; 

> we:=7.2921 15e-5; 

M- := 398618. 
cosi := .5446390348 
h := 59.77 
we := .00007292115 

> K:=evalf(2*Pi/1 1 8*1 e-6*3.3E-2); " 

. K :=. 1757161993 10‘ 8 

> pi :=a A 2*(1 -2*we*sqrt(a*a*a/mu)*cosi); ~ 

pi := a 2 (l -.1258096100 10 ~ 6 Ja*) 

> p2:=(Bessell(0,u)+2*e*Bessell(1 ,u)); ~ ~ “ 

p2 := Bessell( 0, u ) + 2 e Bessell( 1, u ) 

> pe:=exp((6678.1 6-a)/h); ~ " * 

(111.7309687 - .01673080141 a) 

pe := e 

> u:=evalf(simplify(a*e/h)); ~ 

u := .01673080141 a e 

> p:=K*p1*p2*pe; ^ " 

p := .1757161993 10' 8 a 1 (l - .1258096100 10 ~ 6 Ja?) 

(Bessell(0, .01673080141 a e) + 2 e BesseH( 1, .01673080141 a e)) 

(111.7309687 -.01673080141 a) 
e 

> a:=ri)*(1+as*epsilon); — 

> e:=epsilon*es; 

> rO:=6675.16; 

a := rO ( 1 +as e) 
e := 8 es 
rO:= 6675.16 

> p1t:=convert(taylor(p1 , epsilon, 3), poly nom); ~ 

pit := .4150051988 10 8 + .7841517805 10 8 as e + .3118233101 10 8 as 2 e 2 

> epsilon:=0.002740032; ’ ’ 

£ := .002740032 

> pit; 

.4150051988 10 8 + 214860.0971 as + 234.1099365 as 2 
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> pe; 

(.0501924 - .3060089010 os) 
e 

>W, — 

Bessell(0, .00004584293125 (6675.16 + 18.29015201 as) es ) 

+ .005480064 e.s Bessell( 1, .00004584293125 (6675.16 + 18.29015201 as) es) 

> ro'epsiion/n; " 

.3060089010 

> - 
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Appendix I: Reference Tables for Ref 53 


Date: 1 September 1995 

From: Enrico Lorenzini, Smithsonian Astrophysical Observatory 


Reference profiles (last three columns of this file) for 
SEDSAT control law REF53_lSep95 . 

Initial conditions and satellite mass: 

VO - 2.5 m/s, ThetaO - 21 deg (up and fore). Satellite Mass - 36.3 kg 
Reference profiles file: 


Time (s) 

Length (m) 

Speed (m/s) 

0.00 

0.500 

2.500 

8.00 

20.295 

2.492 

16.00 

40.070 

2.485 

24.00 

59.824 

2.478 

32.00 

79.559 

2.472 

40.00 

99.275 

2.466 

48.00 

118.970 

2.461 

56.00 

138.660 

2.456 

64.00 

158.320 

2.451 

72.00 

177.980 

2.447 

80.00 

197.610 

2.443 

88.00 

217.240 

2.440 

96.00 

236.850 

2.437 

104.00 

256.450 

2.434 

112.00 

276.040 

2.432 

120.00 

295.570 

2.430 

128.00 

315.060 

2.429 

136.00 

334.520 

2.428 

144.00 

353.950 

2.427 

152.00 

373.360 

2.427 

160.00 

392.750 

2.428 

168.00 

412.140 

2.428 

176.00 

431.540 

2.430 

184.00 

450.940 

2.431 

192.00 

470.360 

2.434 

200.00 

489.800 

2.436 

208.00 

509.270 

2.439 

216.00 

528.790 

2.443 

224.00 

548.360 

2.447 

232.00 

567.990 

2.451 

240.00 

587.670 

2.455 

248.00 

607.400 

2.461 

256.00 

627.170 

2.466 

264.00 

646.990 

2.472 

272.00 

666.840 

2.479 

280.00 

686.730 

2.486 

288.00 

706.660 

2.493 

296.00 

726.610 

2.501 

304.00 

746.580 

2.509 

312.00 

766.580 

2.518 

320.00 

786.590 

2.527 

328.00 

806.680 

2.537 

336.00 

826.930 

2.547 

344.00 

847.320 

2.557 

352.00 

867.830 

2.568 

360.00 

888.470 

2.578 

368.00 

909.220 

2.590 

376.00 

930.080 

2.601 

384.00 

951.030 

2.613 

392.00 

972.060 

2.626 

400.00 

993.170 

2.639 

408.00 

1014.300 

2.652 


Turns TumRate (1/s) BrakeTums 


3. 3.582 0.000 

32. 3.572 0.000 

60. 3.564 0.000 

88. 3.556 0.000 

117. 3.548 0.000 

1^5. 3.542 0.000 

173. 3.535 0.000 

202. 3.530 0.000 

230. 3.525 0.000 

258. 3.520 0.000 

286. 3.517 0.000 

315. 3.513 0.000 

343. 3.511 0.000 

371. 3.508 0.000 

399. 3.507 0.000 

428. 3.506 0.000 

456. 3.506 0.000 

484. 3.506 0.000 

512. 3.507 0.000 

540. 3.508 0.000 

568. 3.511 0.000 

596. 3.513 0.000 

624. 3.517 0.000 
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0.000 

15313 . 

10.597 

0.000 

15398 . 

10.635 

0.000 

15484 . 

10.675 

0.000 

15569 . 

10.714 

0.000 

15655 . 

10.755 

0.000 

15741 . 

10.795 

0.000 

15828 . 

10.837 

0.000 

15915 . 

10.879 

0.000 

16002 . 

10.921 

0.000 

16089 . 

10.964 

0.000 

16177 . 

11.008 

0.000 

16265 . 

11.052 

0.000 

16354 . 

11.098 

0.000 

16443 . 

11.144 

0.000 

16532 . 

11.190 

0.000 

16623 . 

11.238 

0.000 

16712 . 

11.286 

0.000 

16803 . 

11.335 

0.000 

16894 . 

11.385 

0.000 

16985 . 

11.435 

0.000 

17077 . 

11.486 

0.000 

17168 . 

11.538 

0.000 

17260 . 

11.590 

0.000 

17354 . 

11.644 

0.000 

17447 . 

11.698 

0.000 

17541 . 

11.753 

0.000 

17636 . 

11.809 

0.000 

17731 . 

11.867 

0.000 

17826 . 

11.924 

0.000 

17921 . 

11.983 

0.000 

18017 . 

12.042 

0.000 



2528.00 

10760.000 

5.932 

2536.00 

10807.000 

5.949 

2544.00 

10855.000 

5.966 

2552.00 

10902.000 

5.983 

2560.00 

10950.000 

6.000 

2568.00 

10998.000 

6.018 

2576.00 

11046.000 

6.035 

2584.00 

11095.000 

6.054 

2592.00 

11143.000 

6.072 

2600.00 

11192.000 

6.091 

2608.00 

11241.000 

6.109 

2616.00 

11290.000 

6.129 

2624.00 

11339.000 

6.148 

2632.00 

11388.000 

6.168 

2640.00 

11438.000 

6.188 

2648.00 

11487.000 

6.208 

2656.00 

11537.000 

6.229 

2664.00 

11587.000 

6.249 

2672.00 

11637.000 

6.271 

2680.00 

11687.000 

6.292 

2688.00 

11737.000 

6.313 

2696.00 

11788.000 

6.335 

2704.00 

11839.000 

6.358 

2712.00 

11890.000 

6.380 

2720.00 

11941.000 

6.403 

2728.00 

11992.000 

6.426 

2736.00 

12044.000 

6.449 

2744.00 

12096.000 

6.473 

2752.00 

12148.000 

6.497 

2760.00 

12200.000 

6.521 

2768.00 

12252.000 

6.545 

2776.00 

12304.000 

6.570 

2784.00 

12357.000 

6.595 

2792.00 

12410.000 

6.621 

2800.00 

12463.000 

6.646 

2808.00 

12516.000 

6.672 

2816.00 

12570.000 

6.698 

2824.00 

12623.000 

6.725 

2832.00 

12677.000 

6.751 

2840.00 

12731.000 

6.778 

2848.00 

12786.000 

6.805 

2856.00 

12840.000 

6.833 

2864.00 

12895.000 

6.861 

2872.00 

12950.000 

6.889 

2880.00 

13005.000 

6.917 

2888.00 

13061.000 

6.946 

2896.00 

13116.000 

6.975 

2904.00 

13172.000 

7.004 

2912.00 

13228.000 

7.034 

2920.00 

13285.000 

7.063 

2928.00 

13341.000 

7.093 

2936.00 

13398.000 

7.123 

2944.00 

13455.000 

7.154 

2952.00 

13513.000 

7.185 

2960.00 

13571.000 

7.215 

2968.00 

13628.000 

7.247 

2976.00 

13686.000 

7.278 

2984.00 

13745.000 

7.310 

2992.00 

13803.000 

7.343 

3000.00 

13862.000 

7.375 

3008.00 

13921.000 

7.407 

3016.00 

13980.000 

7.440 

3024.00 

14040.000 

7.473 

3032.00 

14099.000 

7.506 

3040.00 

14159.000 

7.539 

3048.00 

14220.000 

7.573 


18115 . 

12.103 

0.000 

18211 . 

12.165 

0.000 

18309 . 

12.227 

0.000 

18405 . 

12.290 

0.000 

18504 . 

12.354 

0.000 

18603 . 

12.419 

0.000 

18702 . 

12.485 

0.000 

18804 . 

12.553 

0.000 

18903 . 

12.621 

0.000 

19005 . 

12.690 

0.000 

19108 . 

12.761 

0.000 

19210 . 

12.833 

0.000 

19313 . 

12.905 

0.000 

19416 . 

12.979 

0.000 

19521 . 

13.054 

0.000 

19625 . 

13.130 

0.000 

19730 . 

13.207 

0.000 

19837 . 

13.286 

0.000 

19943 . 

13.365 

0.000 

20050 . 

13.446 

0.000 

20157 . 

13.527 

0.000 

20266 . 

13.611 

0.000 

20376 . 

13.696 

0.000 

20486 . 

13.782 

0.000 

20596 . 

13.869 

0.000 

20707 . 

13.957 

0.000 

20820 . 

14.047 

0.000 

20933 . 

14.139 

0.000 

21047 . 

14.232 

0.000 

21161 . 

14.326 

0.000 

21276 . 

14.421 

0.000 

21390 . 

14.518 

0.000 

21508 . 

14.616 

0.000 

21625 . 

14.717 

0.000 

21743 . 

14.818 

0.000 

21862 . 

14.921 

0.000 

21983 . 

15.026 

0.000 

22102 . 

15.132 

0.000 

22223 . 

15.239 

0.000 

22345 . 

15.349 

0.000 

22470 . 

15.461 

0.000 

22593 . 

15.573 

0.000 

22719 . 

15.688 

0.000 

22845 . 

15.805 

0.000 

22971 . 

15.923 

0.000 

23100 . 

16.044 

0.000 

23227 . 

16.165 

0.000 

23357 . 

16.289 

0.000 

23488 . 

16.415 

0.000 

23621 . 

16.543 

0.000 

23752 . 

16.672 

0.000 

23887 . 

16.804 

0.000 

24021 . 

16.937 

0.000 

24159 . 

17.074 

0.000 

24297 . 

17 . 212 

0.000 

24433 . 

17.352 

0.000 

24572 . 

17.494 

0.000 

24715 . 

17.640 

0.000 

24855 . 

17.787 

0.000 

24998 . 

17.937 

0.000 

25142 . 

18.088 

0.000 

25286 . 

18.242 

0.000 

25434 . 

18.398 

0.000 

25579 . 

18.556 

0.000 

25728 . 

18.717 

0.000 

25880 . 

18.882 

0.000 



3056.00 

14281.000 

7.606 

3064.00 

14342.000 

7.641 

3072.00 

14403.000 

7.675 

3080.00 

14465.000 

7.710 

3088.00 

14527.000 

7.745 

3096.00 

14590.000 

7.780 

3104.00 

14652.000 

7.816 

3112.00 

14715.000 

7.852 

3120.00 

14778.000 

7.888 

3128.00 

14842.000 

7.924 

3136.00 

14905.000 

7.960 

3144.00 

14969.000 

7.997 

3152.00 

15033.000 

8.033 

3160.00 

15097.000 

8.069 

3168.00 

15161.000 

8.106 

3176.00 

15226.000 

8.142 

3184.00 

15291.000 

8.178 

3192.00 

15357.000 

8.214 

3200.00 

15423.000 

8.250 

3208.00 

15489.000 

8.285 

3216.00 

15555.000 

8.318 

3224.00 

15622.000 

8.347 

3232.00 

15689.000 

8.373 

3240.00 

15756.000 

8.395 

3248.00 

15823.000 

8.413 

3256.00 

15890.000 

8.426 

3264.00 

15958.000 

8.436 

3272.00 

16025.000 

8.440 

3280.00 

16093.000 

8.440 

3288.00 

16160.000 

8.435 

3296.00 

16228.000 

8.425 

3304.00 

16295.000 

8.409 

3312.00 

16362.000 

8.388 

3320.00 

16429.000 

8.362 

3328.00 

16496.000 

8.329 

3336.00 

16563.000 

8.291 

3344.00 

16629.000 

8.247 

3352.00 

16695.000 

8.197 

3360.00 

16760.000 

8.141 

3368.00 

16825.000 

8.078 

3376.00 

16889.000 

8.010 

3384.00 

16953.000 

7.936 

3392.00 

17016.000 

7.857 

3400.00 

17079.000 

7.771 

3408.00 

17141.000 

7.680 

3416.00 

17202.000 

7.584 

3424.00 

17262.000 

7.483 

3432.00 

17321.000 

7.377 

3440.00 

17380.000 

7.266 

3448.00 

17438.000 

7.152 

3456.00 

17494.000 

7.033 

3464.00 

17550.000 

6.911 

3472.00 

17605.000 

6.786 

3480.00 

17659.000 

6.659 

3488.00 

17711.000 

6.529 

3496.00 

17763.000 

6.397 

3504.00 

17814.000 

6.264 

3512.00 

17863.000 

6.129 

3520.00 

17912.000 

5.994 

3528.00 

17959.000 

5.859 

3536.00 

18006.000 

5.724 

3544.00 

18051.000 

5.589 

3552.00 

18095.000 

5.456 

3560.00 

18138.000 

5.323 

3568.00 

18180.000 

5.192 

3576.00 

18221.000 

5.063 


26032 . 

19.050 

0.000 

26185 . 

19.219 

0.000 

26339 . 

19.392 

0.000 

26496 . 

19.569 

0.000 

26654 . 

19.749 

0.000 

26815 . 

19.933 

0.000 

26974 . 

20.119 

0.000 

27136 . 

20.309 

0.000 

27300 . 

20.502 

0.000 

27466 . 

20.699 

0.000 

27631 . 

20.897 

0.000 

27800 . 

21.099 

0.000 

27969 . 

21.305 

0.000 

28139 . 

21.512 

0.000 

28310 . 

21.723 

0.000 

28485 . 

21.938 

0.000 

28661 . 

22.156 

0.000 

28840 . 

22.378 

0.000 

29020 . 

22.603 

0.000 

29202 . 

22.830 

0.046 

29384 . 

23.053 

0.091 

29570 . 

23.272 

0.136 

29758 . 

23.484 

0.181 

29946 . 

23.689 

0.226 

30136 . 

23.886 

0.271 

30327 . 

24.074 

0.316 

30521 . 

24.256 

0.360 

30715 . 

24.425 

0.404 

30912 . 

24.586 

0.448 

31108 . 

24.733 

0.492 

31308 . 

24.871 

0.535 

31507 . 

24.993 

0.578 

31706 . 

25.101 

0.622 

31908 . 

25.196 

0.664 

32110 . 

25.275 

0.707 

32314 . 

25.339 

0.750 

32517 . 

25.384 

0.792 

32720 . 

25.413 

0.834 

32923 . 

25.422 

0.876 

33126 . 

25.413 

0.918 

33329 . 

25.384 

0.960 

33532 . 

25.337 

1.001 

33734 . 

25.270 

1.042 

33937 . 

25.184 

1.083 

34139 . 

25.078 

1.124 

34339 . 

24.951 

1.165 

34537 . 

24 . 803 

1.205 

34733 . 

24.636 

1.245 

34931 . 

24.452 

1.286 

35127 . 

24.250 

1.325 

35318 . 

24.026 

1.365 

35510 . 

23.789 

1.405 

35700 . 

23.535 

1.444 

35888 . 

23.267 

1.483 

36070 . 

22.980 

1.522 

36254 . 

22.685 

1.561 

36435 . 

22.377 

1.599 

36611 . 

22.056 

1.637 

36788 . 

21.729 

1.676 

36959 . 

21.391 

1.713 

37131 . 

21.050 

1.751 

37297 . 

20.700 

1.789 

37461 . 

20.347 

1.826 

37622 . 

19 . 991 

1.863 

37780 . 

19.633 

1.900 

37936 . 

19.275 

1.937 



3584.00 

18261.000 

4.936 

38088 . 

18.917 

1.974 

3592.00 

18300.000 

4.812 

38238 . 

18.562 

2.010 

3600.00 

18338.000 

4.689 

38385 . 

18.208 

2.046 

3608.00 

18375.000 

4.570 

38529 . 

17.859 

2.082 

3616.00 

18411.000 

4.453 

38671 . 

17.514 

2.118 

3624.00 

18446.000 

4.339 

38809 . 

17.173 

2.154 

3632.00 

18481.000 

4.228 

38948 . 

16.841 

2.189 

3640.00 

18514.000 

4.121 

39079 . 

16.513 

2.224 

3648.00 

18547.000 

4.016 

39212 . 

16.193 

2.260 

3656.00 

18578.000 

3.915 

39337 . 

15.877 

2.294 

3664.00 

18609.000 

3.818 

39464 . 

15.572 

2.329 

3672.00 

18639.000 

3.723 

39586 . 

15.274 

2.364 

3680.00 

18669.000 

3.632 

39710 . 

14.987 

2.398 

3688.00 

18698.000 

3.544 

39830 . 

14.707 

2.432 

3696.00 

18726.000 

3.459 

39946 . 

14.434 

2.466 

3704.00 

18753.000 

3.377 

40059 . 

14.169 

2.500 

3712.00 

18780.000 

3.298 

40173 . 

13.914 

2.533 

3720.00 

18806.000 

3.222 

40283 . 

13.667 

2.566 

3728.00 

18831.000 

3.149 

40389 . 

13.427 

2.600 

3736.00 

18856.000 

3.079 

40496 . 

13.196 

2.633 

3744.00 

18880.000 

3.011 

40599 . 

12.973 

2.665 

3752.00 

18904.000 

2.946 

40703 . 

12.758 

2.698 

3760.00 

18928.000 

2.884 

40807 . 

12.553 

2.730 

3768.00 

18950.000 

2.824 

40903 . 

12.351 

2.762 

3776.00 

18973.000 

2.766 

41004 . 

12.161 

2.793 

3784.00 

18995.000 

2.711 

41101 . 

11.976 

2.825 

3792.00 

19016.000 

2.658 

41194 . 

11.797 

2.856 

3800.00 

19037.000 

2.607 

41287 . 

11.626 

2.886 

3808.00 

19058.000 

2.558 

41381 . 

11.463 

2.917 

3816.00 

19078.000 

2.511 

41471 . 

11.304 

2.947 

3824.00 

19098.000 

2.466 

41561 . 

11.153 

2.977 

3832.00 

19118.000 

2.422 

41652 . 

11.008 

3.007 

3840.00 

19137.000 

2.380 

41738 . 

10.867 

3.036 

3848.00 

19156.000 

2.340 

41825 . 

10.731 

3.066 

3856.00 

19174.000 

2.301 

41908 . 

10.600 

3.095 

3864.00 

19193.000 

2.264 

41996 . 

10.476 

3.123 

3872.00 

19210.000 

2.228 

42075 . 

10.352 

3.152 

3880.00 

19228.000 

2.193 

42158 . 

10.235 

3.180 

3888.00 

19246.000 

2.159 

42243 . 

10.124 

3.208 

3896.00 

19263.000 

2.126 

42323 . 

10.013 

3.236 

3904.00 

19280.000 

2.095 

42403 . 

9.908 

3.264 

3912.00 

19296.000 

2.064 

42479 . 

9.803 

3.291 

3920.00 

19313.000 

2.034 

42560 . 

9.704 

3.318 

3928.00 

19329.000 

2.005 

42636 . 

9.606 

3.346 

3936.00 

19345.000 

1.977 

42713 . 

9.511 

3.372 

3944.00 

19360.000 

1.950 

42785 . 

9.417 

3.399 

3952.00 

19376.000 

1.923 

42863 . 

9.328 

3.426 

3960.00 

19391.000 

1.897 

42936 . 

9.239 

3.452 

3968.00 

19406.000 

1.872 

43009 . 

9.153 

3.478 

3976.00 

19421.000 

1.847 

43082 . 

9.069 

3.504 

3984.00 

19436.000 

1.822 

43156 . 

8.986 

3.530 

3992.00 

19450.000 

1.798 

43225 . 

8.904 

3.556 

4000.00 

19465.000 

1.775 

43300 . 

8.825 

3.581 

4008.00 

19479.000 

1.752 

43369 . 

8.746 

3.607 

4016.00 

19493.000 

1.730 

43439 . 

8.668 

3.632 

4024.00 

19506.000 

1.707 

43505 . 

8.589 

3.657 

4032.00 

19520.000 

1.686 

43575 . 

8.515 

3.682 

4040.00 

19533.000 

1.664 

43641 . 

8.438 

3.707 

4048.00 

19547.000 

1.643 

43712 . 

8.365 

3.732 

4056.00 

19560.000 

1.622 

43779 . 

8.291 

3.756 

4064.00 

19573.000 

1.601 

43845 . 

8.217 

3.781 

4072.00 

19585.000 

1.581 

43907 . 

8.142 

3.805 

4080.00 

19598.000 

1.561 

43974 . 

8.070 

3.830 

4088.00 

19610.000 

1.541 

44036 . 

7.995 

3.854 

4096.00 

19622.000 

1.521 

44098 . 

7.921 

3.878 

4104.00 

19635.000 

1.501 

44166 . 

7.851 

3.902 



4112.00 

19646.000 

1.481 

44224 . 

7.774 

3 # 926 

4120.00 

19658.000 

1.462 

44287 . 

7.701 

3.950 

4128.00 

19670.000 

1.442 

44350 . 

7.628 

3.974 

4136.00 

19681.000 

1.423 

44409 . 

7.551 

3.998 

4144.00 

19693.000 

1.403 

44472 . 

7.477 

4.021 

4152.00 

19704.000 

1.384 

44531 . 

7.401 

4.045 

4160.00 

19715.000 

1.365 

44590 . 

7.323 

4.069 

4168.00 

19726.000 

1.345 

44649 . 

7.246 

4.092 

4176.00 

19736.000 

1.326 

44703 . 

7.165 

4.116 

4184.00 

19747.000 

1.306 

44763 . 

7.085 

4.139 

4192.00 

19757.000 

1.286 

44817 . 

7.003 

4.162 

4200.00 

19767.000 

1.267 

44872 . 

6.918 

4.186 

4208.00 

19777.000 

1.247 

44926 . 

6.833 

4 . 209 

4216.00 

19787.000 

1.227 

44981 . 

6.747 

4.233 

4224.00 

19797.000 

1.207 

45036 . 

6.659 

4 . 256 

4232.00 

19807.000 

1.186 

45092 . 

6.569 

4 . 279 

4240.00 

19816.000 

1.166 

45142 . 

6.475 

4.303 

4248.00 

19825.000 

1.145 

45192 . 

6.380 

4.326 

4256.00 

19834.000 

1.123 

45242 . 

6.282 

4.349 

4264.00 

19843.000 

1.102 

45292 . 

6.181 

4.373 

4272.00 

19852.000 

1.080 

45343 . 

6.079 

4.396 

4280.00 

19861.000 

1.058 

45394 . 

5.973 

4 .420 

4288.00 

19869.000 

1.035 

45439 . 

5.862 

4.443 

4296.00 

19877.000 

1.012 

45484 . 

5.748 

4.467 

4304.00 

19885.000 

0.988 

45530 . 

5.628 

4.491 

4312.00 

19893.000 

0.963 

45575 . 

5.503 

4 . 515 

4320.00 

19901.000 

0.938 

45621 . 

5.373 

4.539 

4328.00 

19908.000 

0.911 

45661 . 

5.235 

4.564 

4336.00 

19915.000 

0.884 

45702 . 

5.092 

4.588 

4344.00 

19922.000 

0.856 

45742 . 

4.942 

4 . 613 

4352.00 

19929.000 

0.826 

45782 . 

4.787 

4.637 

4360.00 

19935.000 

0.796 

45817 . 

4.624 

4.662 

4368.00 

19942.000 

0.765 

45858 . 

4.456 
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